function demod = asynchdemod(r)
%Takes a signal r and demodulates using asynchronous modulation. 
%Assumptions:
% The input signal is sampled at 88.2kHz.
% The input signal is 10 seconds long.
% The carrier frequency is 15kHz. 
%Returns:
% The demodulated signal downsampled to 44.1kHz.

%Overall sampling rate, 88.2kHz.
rate = 44100*2;

%Filter order
order = 50;

%A lowpass filter, operates at a sampling rate of 88.2kHz, with cutoff 8kHz.
b = fir1(order, 8e3/(rate/2)); a=1;

%Time axis, our signals are exactly 10 seconds long
t = 0:1/rate:10;
t = t(1:length(t)-1)';

%Plot spectrum of r
%We need a high fft_len to get good resolution
fft_len = 2^16;
R = fft(r,fft_len);
mag_R = abs(R)/fft_len;

%Our frequency scale, in Hz, positive frequencies only
f = rate*(0:fft_len/2)/fft_len;

figure(1);
plot(f,mag_R(1:fft_len/2+1))
title('Spectrum of AM signal');
xlabel('Frequency (Hz)');
ylabel('Magnitude');

%Asynchronous demodulation
%Multiply by sin/cos
c_15k = cos(2*pi*15000*t);
s_15k = sin(2*pi*15000*t);

xc = c_15k.*r;
xs = s_15k.*r;

% LPF at carrier frequency
lpfc = fir1(order, 15000/(rate/2));

xc = filter(lpfc,a,xc);
xs = filter(lpfc,a,xs);

% Throw out first (order) samples of xc and yc because the filter
% is not effective for those points

xc = xc(order:end);
xs = xs(order:end);

% Square, sum, and square root

s= ((xc.^2)+(xs.^2)).^(.5);

% Subtract average to remove constant added by injection of carrier
s = s-mean(s);

% LPF
demod = filter(b,a, s);
%Rescale
demod = 2*demod;
%Downsample to the audio range
demod = decimate(demod,2);

%demod is now sampled at 44.1kHz
%Play sound 
soundsc2(demod,44100);
