Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % 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');
Add Comment
Please, Sign In to add comment