% acquire signal
[s,fs,nbits] = wavread('rexkwondo.wav');
% low-pass filter
f_cutoff = 1e3/fs;
W_trans = 100/fs;
b = firpm(51, 2*[0 f_cutoff-W_trans f_cutoff+W_trans 0.5], [1 1 0 0]);
sn = filter(b,[1],s);
sn=sn-mean(sn);  % remove DC component.

% define time and frequency axes
tt = (0:length(sn)-1)/fs;
tt = tt';
ff = (0:length(sn)-1)*fs/length(sn) - fs/2;

%plot original signal and spectrum
figure(1)

subplot(2,1,1)
plot(tt,sn)
xlabel('Time(s)')
ylabel('Amplitude')

Sn = fft(sn);
subplot(2,1,2)
plot(ff,abs(fftshift(Sn)));
xlabel('Frequency (Hz)')
ylabel('|Fourier Transform|')

p=audioplayer(sn,fs);
play(p)

pause

% Multiplication with complex exponential

fc = 1500;
x1 = sn.*exp(j*2*pi*fc*tt);

%plot  signal and spectrum
figure(2)

subplot(2,1,1)
plot(tt,x1)
xlabel('Time(s)')
ylabel('Amplitude')

X1 = fft(x1);
subplot(2,1,2)
plot(ff,abs(fftshift(X1)));
xlabel('Frequency (Hz)')
ylabel('|Fourier Transform|')

p1=audioplayer(x1,fs);
play(p1)

pause

% Multiplication with sinusoid 
x = sn.*cos(2*pi*fc*tt);

%plot signal and spectrum
figure(3)

subplot(2,1,1)
plot(tt,x)
xlabel('Time(s)')
ylabel('Amplitude')

X = fft(x);
subplot(2,1,2)
plot(ff,abs(fftshift(X)));
xlabel('Frequency (Hz)')
ylabel('|Fourier Transform|')

px=audioplayer(x,fs);
play(px)

pause

% Demod: multiply again
y = x.*cos(2*pi*fc*tt);

%plot signal and spectrum
figure(4)

subplot(2,1,1)
plot(tt,y)
xlabel('Time(s)')
ylabel('Amplitude')

Y = fft(y);
subplot(2,1,2)
plot(ff,abs(fftshift(Y)));
xlabel('Frequency (Hz)')
ylabel('|Fourier Transform|')

py=audioplayer(y,fs);
play(py)

pause

% lowpass filter
yn = filter(b,[1],y);
%plot signal and spectrum
figure(5)

subplot(2,1,1)
plot(tt,yn)
xlabel('Time(s)')
ylabel('Amplitude')

Yn = fft(yn);
subplot(2,1,2)
plot(ff,abs(fftshift(Yn)));
xlabel('Frequency (Hz)')
ylabel('|Fourier Transform|')

pyn=audioplayer(yn,fs);
play(pyn)

pause

