Saturday, July 18, 2009

Finite impulse response (FIR) filter design

% % prog7b1.m, p442 Program name
%
fc=0.53; % cutoff frequency (normalized to Fs/2)
N= 10; % Filter length (number of taps)
hd=fir1(N-1,fc,boxcar(N)); % Calculate truncated ideal impulse response coeffs
wn=hamming(N); % Calculate Hamming window coefficients
hn=fir1(N-1,fc,wn); % Obtain windowed coefficients

===============================
% % prog7b2.m, p444 File name %
FS=1000; % Sampling frequency
FN=FS/2; % Nyquist frequency
N=73; % Filter length
beta=5.65; % Kaiser ripple parameter

fc1=125/FN; % Normalized cut off frequencies
fc2=275/FN;
FC=[fc1 fc2]; % Band edge frequencies
hn=fir1(N-1, FC, kaiser(N, beta)); %Obtain windowed filter coefficients
[H,f]=freqz(hn, 1, 512, FS); % Compute frequency response
mag=20*log10(abs(H));
plot(f, mag), grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude Response (dB)')

===============================
% % Program name - prog7b3.m, p446
% m-file for calculating optimal FIR filter coefficients and
% plotting frequency response

%
Fs=10000; % Sampling frequency
N=41; % Filter length
WT=[10 3 10]; % Weights of the deviations in the bands
Hd=[0 0 1 1 0 0]; % Desired magnitude response in the bands
F=[0 0.1 0.2 0.3 0.4 1]; % Band edge frequencies
b = remez(N-1, F, Hd, WT); % Compute the filter coefficients
[H, f] = freqz(b, 1, 512, Fs);% Compute the frequency response
mag = 20*log10(abs(H)); % of filter and plot it
plot(f, mag), grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')

=================================
% % Program name - prog7b4a.m, p449.
% m-file for calculating the optimal FIR filter coefficients
% and plotting frequency response for Example 7B.4

%
Fs=50000; % Sampling frequency
% Filter length
Ap=1; % Pass band ripple in dB

As=45; % Stop band attenuation in dB
M=[0 1 0]; % Desired magnitude response
F=[10000, 12000, 16000, 18000] ; % Band edge frequencies
dp=(10^(Ap/20)-1)/(10^(Ap/20)+1); % Pass and stop band ripples
ds=10^(-As/20);
dev=[ds dp ds];
[N1, F0, M0, W] = remezord(F, M, dev, Fs);% Determine filter order
[b delta] = remez(N1, F0, M0, W);% Compute the filter coefficients
[H, f] = freqz(b, 1, 1024, Fs); % Compute the frequency response
mag = 20*log10(abs(H)); % of filter and plot it
plot(f, mag), grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')

===========================================
% % Program name - prog7b4b.m, p450.
% An alternative m-file for calculating the optimal FIR filter
% coefficients and plotting frequency response for Example 7B.4

%
N=44;
Fs=50000; % Sampling frequency
% Filter length

Ap=1; % Pass band ripple in dB
As=45; % Stop band attenuation in dB
M=[0 0 1 1 0 0]; % Desired magnitude response
F=[0, 0.4, 0.48, 0.64, 0.72 1] ; % Band edge frequencies
dp=(10^(Ap/20)-1)/(10^(Ap/20)+1);
ds=10^(-As/20);
W=[dp/ds, 1, dp/ds];
dev=[ds ds dp dp ds ds ];
[b delta] = remez(N-1, F, M, W);% Compute the filter coefficients
[H, f] = freqz(b, 1, 1024, Fs); % Compute the frequency response
mag = 20*log10(abs(H)); % of filter and plot it
plot(f, mag), grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')

==============================================
% % Program name - prog7b5.m, p451.
% m-file for computing the coefficients of an FIR frequency sampling

% filter
%
Fs=2000; % Sampling frequency
N=15; % Filter length
fd=[0 1/7 2/7 3/7 4/7 5/7 6/7 1]; % Frequency sampling points
Hd=[1 1 1 1 0.5571 0.0841 0 0]; % Frequency samples
hn=fir2(N-1, fd, Hd); % Compute the impulse response coeffs.
[H, f] = freqz(hn, 1, 512, Fs); % Plot the magnitude frequency response
plot(f, abs(H)), grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude ')

==================================================
% % Program name - prog7b6.m, p453.
% m-file for calculating the coefficients of an FIR filter
% with arbitrary magnitude response
%
Fs=2000;
N=110;
fd=[0 0.15 0.25 0.45 0.5 0.75 0.85 1]; % Frequency sampling points
Hd=[1 1 0.3 0.3 0.1 0.1 0 0]; % Frequency samples
hn=fir2(N-1, fd, Hd); % Compute the impulse response
[H, f] = freqz(hn, 1, 512, Fs); % Plot the magnitude response
plot(f, abs(H)), grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude ')

=================================================
function Filt_win
%This program uses the window method to design FIR filters.
%The following windows are supported: Rectangular, Hamming, Hanning, Blackman,
%and Kaiser.

% ftype: [] for lowpass filter
% 'high' for highpass filter
% [] for bandpass filter
% 'stop' for bandstop filter
% wtype: 1 for Rectangular window
% 2 for Hamming window
% 3 for Hanning window
% 4 for Blackman window
% 5 for Kaiser window

clear all;
wtype = 2 ;
Fs = 8000; %sampling frequency
n = 52; %the order of the filter
Wn = 0.4375; %normalised cutoff frequency (2*Fstop/Fs)
beta = 0; %Kaiser window beta parameter
ftype=[];
window=[]; %Hamming window will be used if variable window is empty

if wtype==1
window=ones(1,n+1); %Rectangular window
elseif wtype==3
window=hanning(n+1); %Hanning window
elseif wtype==4
window=blackman(n+1); %Blackman window
elseif wtype==5
window=kaiser(n+1,beta);%Kaiser window
end
b = fir1(n,Wn,ftype,window)

% Plot the frequency response
freqz(b,1,1024,Fs);

=========================================
function Firfilt
%FIR filtering program - program performs FIR filtering using coefficients
%and data in user specified files.
clear all;
data_file = fopen('filterin.dat','r'); %read the data to be filtered from file

x = fscanf(data_file,'%f');
coef_file = fopen('filtcoef.dat','r'); %read the filter coefficients from file
[b,n]=fscanf(coef_file,'%f',inf);
fclose('all');

y=filter(b,1,x); %filtering

%Plot the results
bx=fft(x,512);
by=fft(y,512);
[h, w]=freqz(b,1,1024);
subplot(3,1,1),plot(1000*(0:255)/512,abs(bx(1:256))), ylabel('X(f)');
subplot(3,1,2),plot(w/pi,abs(h)), ylabel('H(f)');
subplot(3,1,3),plot(1000*(0:255)/512,abs(by(1:256))), ylabel('Y(f)');

===================================================================
function Fresamp
%Design arbitrary FIR filter by frequency sampling method
n=15;
m = [1 1 1 1 0.5571 0.0841 0 0];
f = [0 1/7 2/7 3/7 4/7 5/7 6/7 1];
window = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; % n+1 elements long

b = fir2(n,f,m,window); % get filter coefficients

%Output the results
disp('The coefficients are:'); disp(b);
[h,w] = freqz(b,1,128);
plot(f,m,w/pi,abs(h)); title('Desired and Actual Frequency Responses');


==================================================
% % Program name - REMEZ1.m
%
Fs=15000; % sampling frequency

N=41; % Filter length
WT=[10 3 10]; % Weights of the deviations in the bands
Hd=[0 0 1 1 0 0]; % Desired magnitude response in the bands
F=[0 0.06 0.12 0.14667 0.206667 1]; % Band edge frequencies
[b, delta] = remez(N-1, F, Hd, WT); % Compute the filter coefficients
[H, f] = freqz(b, 1, 512, Fs); % Compute the frequency response
mag = 20*log10(abs(H)); % of filter and plot it
plot(f, mag)
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')

================================================
% % Program name – REMEZex2.m
%
N=44
Fs=50000; % Sampling frequency
% Filter length

Ap=1; % Pass band ripple in dB
As=45; % Stop band attenuation in dB
M=[0 0 1 1 0 0]; % Desired magnitude response
F=[0, 0.4, 0.48, 0.64, 0.72 1] ; % Band edge frequencies
dp=(10^(Ap/20)-1)/(10^(Ap/20)+1);
ds=10^(-As/20);
W=[dp/ds, 1, dp/ds];
dev=[ds ds dp dp ds ds ];
[b delta] = remez(N-1, F, M, W);% Compute the filter coefficients
[H, f] = freqz(b, 1, 1024, Fs); % Compute the frequency response
mag = 20*log10(abs(H)); % of filter and plot it
plot(f, mag), grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
================================================

No comments: