function[N] = bellangord(dp, ds, Fp, Fs, FT)
num = -2*log10(10*dp*ds);
den = 3*(Fs - Fp)/FT;
N = ceil(num/den)-1;
function[N] = kaiord(dp, ds, Fp, Fs, FT);
num = -20*log10(sqrt(dp*ds))-13;
den = 14.6*(Fs - Fp)/FT;
N = ceil(num/den);
function [y,ssp,iter] = minphase(g);
% Extracts the minimum phase factor y from a linear-phase filter.
% Input: g = [g(0) g(1) . . . g(N)] (row vector) where
% the g vector is the right-half of a linear-phase FIR filter.
% Assumption: Unit-circle zeros of g are of even multiplicity.
% Copyright (c) January 2002 by H.J. Orchard and A.N. Willson, Jr.
% Initialize y (poly. with all zeros at z = 0)
y = [1 zeros(1,length(g)-1)];
ssp = realmax;% A large number (for previous norm)
ss = ssp/2; % Smaller large no. (for current norm)
iter = 0; e = 0; % Initialize iter. counter & correction vector
while ( ss <>
y = y + e'; ssp = ss;% Update y and move old norm(e) value
iter = iter + 1; % Increment the iteration count
AR = toeplitz([y(1),zeros(1, length(h)-1)], y);
M = toeplitz([y(length(g)),zeros(1,length(g)-1)],fliplr(y))
AL = fliplr(M);
A = AL + AR;% Create the A matrix
b = g' - AL*y'; % and the b vector
e = A\b;% Solve Ae = b for the correction vector e
ss = norm(e); % Get norm to see if it is still decreasing
end;
% Program 10_1
% Estimation of FIR Filter Order Using remezord
%
fedge = input('Type in the bandedges = ');
mval = input('Desired magnitude values in each band = ');
dev = input('Allowable deviation in each band = ');
FT = input('Type in the sampling frequency = ');
[N, fpts, mag, wt] = remezord(fedge, mval, dev, FT);
fprintf('Filter order is %d \n',N);
% Program 10_2
% Design of Equiripple Linear-Phase FIR Filters
%
format long
fedge = input('Band edges in Hz = ');
mval = input('Desired magnitude values in each band = ');
dev = input('Desired ripple in each band =');
FT = input('Sampling frequency in Hz = ');
[N,fpts,mag,wt] = remezord(fedge,mval,dev,FT);
b = remez(N,fpts,mag,wt);
disp('FIR Filter Coefficients'); disp(b)
[h,w] = freqz(b,1,256);
plot(w/pi,20*log10(abs(h)));grid;
xlabel('\omega/\pi'); ylabel('Gain, dB');
% Program 10_3
% Design of a minimum-phase lowpass FIR filter
%
Wp = 0.45; Ws = 0.6; Rp = 2; Rs = 26;
% Desired ripple values of minimum-phase filter
dp = 1- 10^(-Rp/20); ds = 10^(-Rs/20);
% Compute ripple values of prototype linear-phase filter
Ds = (ds*ds)/(2 - ds*ds);
Dp = (1 + Ds)*((dp + 1)*(dp + 1) - 1);
% Estimate filter order
[N,fpts,mag,wt] = remezord([Wp Ws], [1 0], [Dp Ds]);
% Design the prototype linear-phase filter H(z)
[b,err,res] = remez(N, fpts, mag, wt);
K = N/2;
b1 = b(1:K);
% Design the linear-phase filter G(z)
c = [b1 (b(K+1) + res.error(length(res.error))) fliplr(b1)]/(1+Ds);
zplane(c) % Plot the zeros of G(z)
pause
c1 = c(K+1:N+1);
[y, ssp, iter] = minphase(c1);
zplane(y) % Plot the zeros of the minimum-phase filter
pause
[hh,w] = freqz(y,1, 512);
% Plot the gain response of the minimum-phase filter
plot(w/pi, 20*log10(abs(hh)));grid
xlabel('\omega/\pi');ylabel('Gain, dB');
% Program 10_4
% Kaiser Window Generation
%
fpts = input('Type in the bandedges = ');
mag = input('Type in the desired magnitude values = ');
dev = input('Type in the ripples in each band = ');
[N,Wn,beta,ftype] = kaiserord(fpts,mag,dev)
w = kaiser(N+1,beta); w = w/sum(w);
[h,omega] = freqz(w,1,256);
plot(omega/pi,20*log10(abs(h)));grid;
xlabel('\omega/\pi'); ylabel('Gain, dB');
% Program 10_5
% Lowpass Filter Design Using the Kaiser Window
%
fpts = input('Type in the bandedges = ');
mag = input('Type in the desired magnitude values = ');
dev = input('Type in the ripples in each band = ');
[N,Wn,beta,ftype] = kaiserord(fpts,mag,dev)
kw = kaiser(N+1,beta);
b = fir1(N,Wn, kw);
[h,omega] = freqz(b,1,512);
plot(omega/pi,20*log10(abs(h)));grid;
xlabel('\omega/\pi'); ylabel('Gain, dB');
% Program 10_6
% Design of Multiband FIR Filter Using Hamming Window
%
fpts = [0 0.28 0.3 0.5 0.52 1];
mval = [0.3 0.3 1.0 1.0 0.7 0.7];
b = fir2(100,fpts,mval);
[h,omega] = freqz(b,1,512);
plot(omega/pi,abs(h));grid;
xlabel('\omega/\pi'); ylabel('Magnitude');
No comments:
Post a Comment