% from demo program filtdem in lab1
%% create the test signal
clear all;
Fs = 200; % sampling frequency
t=(1:100)/Fs ; % obtain sampling time vector [ 1*1/Ts  2*1/Ts ...  100/Ts]
s1= sin(2*pi*10*t);   % create 5 hz test signal , 100 samples
s2= sin(2*pi*30*t);%s2= sin(2*pi*15*t);   % create 15 hz test signal , 100 samples
s3= sin(2*pi*50*t);%s3= sin(2*pi*30*t);   % create 30 hz test signal , 100 samples

s= s1+s2+s3 ; % Create a test signal that contain 3 different frequecy 5,15 and 30Hz
plot(t,s)  % plot the signal vs time(seconds)

%% Design the filter

%b = fir1( n , wn );
a=-1;
b=fir1(33,0.2);
%figure()
%freqz(b,a,512,Fs);
% implement an IIR bandpass filter using the elliptic method
%[b,a] = ellip(4,0.1,40,[10 20]*2/Fs); 
%From ellip help instruction, cut off frequency Wp must be 0.0 < Wp < 1.0, 
% with 1.0 corresponding to half the sample rate, thus frequency used in this function is divided by Fs/2

% check the frequency response of the filter that has been designed earlier
% to see if it meet the desired specification
figure()
freqz(b,a,512,Fs); % short cut way to get frequency response plot (Fs = sampling frequency in Hz)
figure();
[H,w] = freqz(b,a,512);
%  Obtain the frequency response H, evaluated at N=512 points equally spaced around the
%  upper half of the unit circle, frequency ranging from 0 to pi
% w in radians/sample (discrete time frequency)
f = w*Fs/(2*pi);  %conversion to frequency in Hz based on formula w = 2*pi*f/Fs

% Plot frequency response in Hz
figure,plot(f,abs(H)), title('frequency response of the elliptic IIR bandpass filter');
ylabel('magnitude of H(w)'); xlabel('Frequency in Hz'); 

% apply the filter to filter the input signal s and get output signal sf
sf = filter(b,a,s);
figure, plot(t,sf);
xlabel('Time in seconds');
ylabel('waveform magnitude');
axis([0 1 -1 1]);


%% Find the frequency content (freq. spectrum) of the signal before filtering

% Obtain freq spectrum by performing 512 point DFT using fft algorithm, a higher number gives spectrum
% with better resolution, 
% Take note that the DFT sample below are for frequency ranging from 0->fsampling/2
S = fft(s,512) ;   % DFT samples vector [S(1) S(2) .....S(512) ]

% Find freq spectrum for signal after filtering
SF = fft(sf,512) ; % DFT samples vector [SF(1) SF(2) .....SF(512) ]


% Plot the frequency spectrum S for the signal before filtering  
% frequency in hertz
% Due to the periodic nature of the frequency spectrum of discrete time signal
% we only take half of the 512 samples
wk = (0:511) * (2*pi/512)  ;  % frequency vector in omega 0-2pi
fk = wk * Fs/(2*pi);
figure, plot(fk(1:256), abs(S(1:256)), 'b-' ); grid on  % black line 
xlabel(' Frequency (Hz)  ');
ylabel('Magnitude of fourier transform samples');

hold on
% Plot the frequency spectrum SF for the filtered signal  
plot(fk(1:256), abs(SF(1:256)), 'r--'); % red dotted line
xlabel(' Frequency (Hz)  ');
ylabel('Magnitude of fourier transform samples');