%%
close all; clear all;
useZoom = 0;

[m, fs] = wavread('audio48kHz.wav');
%sound(m, fs);
Ts = 1/fs;
plotspec(m,Ts)
if(useZoom)
    subplot(2,1,1); xlim([1 1.005])
end
%% Modulation
figure
t = (0:Ts:((length(m)-1)*Ts))';
fc = 10000; 
x = m.*(sqrt(2)*(cos(2*pi*fc*t))); % Direct Multiplication Modulator
%x = m.*(sqrt(2)*(cos(2*pi*fc*t)>=0)); % Switching Modulator
plotspec(x,Ts)
if(useZoom)
    subplot(2,1,1); xlim([1 1.005]); hold on; 
    plot(t,sqrt(2)*m,'g--'); plot(t,-sqrt(2)*m,'r--');
end
%sound(x, fs);

%% Demodulation Step 1
figure
v = x.*(sqrt(2)*cos(2*pi*fc*t)); % Direct Multiplication Demodulator
%v = x.*(sqrt(2)*(cos(2*pi*fc*t)>=0)); % Switching Demodulator
plotspec(v,Ts)
%sound(v, fs);
if(useZoom)
    subplot(2,1,1); xlim([1 1.005])
end
%% Demodulation Step 2
figure
% Build the LPF
fcutoff = 0.5e4;
freqs=[0 2*fcutoff/fs 2.2*fcutoff/fs 1];            % Scaled by fs/2
amps=[1 1 0 0];
b=firpm(100,freqs,amps);         % specify the LP filter
mhat=filter(b,1,v);               % do the filtering
plotspec(mhat,Ts)       % plot the output spectrum
%sound(mhat, fs);
if(useZoom)
    subplot(2,1,1); xlim([1 1.005])
    % Note that there is some delay in mhat from the filtering operation
end

