For the ideal H'_d(w) and h_d(n):

where H_d(e^jw) is a periodic function with period = 2p.

The inverse is:

Previous chapter's work suggests that the impulse responses for ideal filters are infinite in length.

So in actuality, we design finite filters by truncation of the form:

If the original h_d(n) had an infinite impulse response,

we have the operation: h(n) = h_d(n) w_R(n)

where w_R(n) = 1 for 0<n<M and = 0 otherwise

So, if we know the original nature of the answer, since multiplication in the time domain is like convolution in the frequency domain:

where:

which is our "old friend"

which we have examined before. This is the convolution that gives us many bumps as it passes by discontinuities in the frequency domain.

This is called the Gibbs phenomenon or boxcar effect.

We therefore try to select windows with preferable frequency responses, i.e.--better represent a delta function or have smaller side lobes.

The windows look just like those we examined in Chapter 7. Only now they are used for a different purpose. With regard to spectrum analysis, the windows were multiplied times the data points x(n), e.g.

~x(n) = x(n) w(n) where w(n) is a window defined for 0 < n < N-1

;

In FIR FILTER DESIGN, the windows used are basically the same, however, they are used differently----they are used to MODIFY THE IMPULSE RESPONSE OF THE SYSTEM, defined for n = 0, 1, ...M rather than to modify the data points (n = 0, 1, N-1). Also, the desirability of each is different with respect to the different applications. For example, in Chapter 7, we indicated that nothing was much better than the HANNING because it had such a sharp roll-off. For designing filters, the height of the largest side-lobe is a more important feature.

Common ones using impulse response ranges are:

Main Lobes for Hanning and Hamming: 8p/M radians wide;

Main Lobes for Blackman: 12p/M radians wide;

FIRS and LINEAR PHASE DELAY --- THE REAL ADVANTAGE

Examples:	h(M-n) = h (n)	for n=0,1,...M
or	h(M-n) = -h (n)	for n=0,1,...M

We are going to focus on CAUSAL, SYMMETRIC TYPE SYSTEMS:

GENERAL EQUATIONS SHIFT TO:

i.e., assumes the impulse response is finite and = 0 outside the 0 to M range.

Two equations we have used have the desired symmetric properties:

ideal lowpass filter:

ideal lowpass differentiator:

ASIDE ON IDEAL DIFFERENTIATION

(Integration by parts Proof on Page 152)

Therefore, if we wish to preserve symmetry, i.e. assuming we begin with a function H_d(ejw) that has the desired even properties: h(n) = h(-n),

and we want to window it with our even windows and have a window of M+1 points:

we first shift:

Then: h_d(n), 0<n<M, will have the desired properties: h_d(n) = h_d(M-n)

ideal lowpass filter:

ideal lowpass differentiator:

DESIGN BY WINDOWING:

RECTANGULAR WINDOW APPLIED TO IDEAL LPF

HANNING WINDOW APPLIED TO IDEAL LP FILTER

1. The only parameter to select for the windows is M.

2. M does not control the height of the sidelobes, only the width of the mainlobe and transition band.

3. d₁ = d₂ in the final filter product.

COMPARISONS TO RECTANGULAR WINDOWS

The Peak to Sidelobe Level is relative to the window itself, i.e. the difference in dBs from W'(0) to the highest sidelobe.

The table value for Transition Bandwidth is the value of "k" to use in the equation given. For example, if M=10 and the window is Hamming, the Transistion Bandwidth is: 2p(3.3)/10=.66p radians

This table therefore gives us the d₁ and d₂ values for each window, therefore d₁, which encompasses the passband is no longer associated with a 3 dB down level for w_c. The w_d level specifies the frequency where the amplitude is at .5 or -3dBs, independent of the window selected and the value of M.

e.g. for Hanning Window: 20 log10(d₁) = -44 implies that
d₁ = 10^(-44/20) = .063

KAISER WINDOWS: This is a family of windows that serves to mimic all of the other windows, plus gives us control of d₂:

It has the general formula:

where: I_o(b) is a zero-order modified Bessel function of the first kind.

Your text gives you a power series formula--
but it is also a MATLAB FUNCTION.

b = 0 is rectangular window, b = 3.86 is comparable to Hanning, etc

FIR Filter design assumes a common d:

i.e. | H(e^jw) | > 1 - d in passband, | H(e^jw) |< d in stopband

for lowpass:

	A = -20log10d2
	= 0.1102(A-8.7)	A>50
ß	= 0.5842(A-21)0.4+.07886 (A-1)	21<A<50
	= 0	A<21
M	=

DESIGN PROBLEM

wc = 0.2p,

wr = 0.3p,

A = 50 dB

With an A of 50, we see that the Hamming (-53 dB sidelobe) will work. Now determine M:

w_r - w_c = 0.1p = 2 p (3.3)/M or M = 2 p (3.3)/ 0.1p = 20 (3.3) = 66

We would select the midpoint between w_c and w_r for w_d

i.e. w_d = 0.25p

Using Kaiser Window:

Compare this to the IIR Butterworth Design where N=15, i.e.

So even using the worst programming practices for the zeros, we would get by with 30 multiplications!

IDEAL DIFFERENTIATER with RECTANGULAR WINDOW

IDEAL DIFFERENTIATER with HANNING WINDOW

NOTE: TO SPEED THINGS UP, I AM SHIFTING TO MATLAB FOURIER TRANSFORM FOR THIS DEMONSTRATION: DATA IS ASSUMED TO START AT n=0---which is now true.

function [X,w] = clwfour(x,nw);

Y=fft(x,nw);

w=[-1:2/nw:1-2/nw];

X=fftshift(Y);

clg;subplot(221);samp(x);xlabel('n');title('real x[n]');...

subplot(223);plot(w,20*log10(abs(X)));xlabel('norm w');title('20*log');...

subplot(222);plot(w,abs(X));v=axis;v(3)=0;axis(v);...

xlabel('norm w');title ('modulus');grid;...

subplot(224);plot(w,angle(X));v=axis;v(3)=-4;v(4)=4;axis(v);...

xlabel('norm w');title('phase');grid;

subplot(111);

CLFIR DEMONSTRATION

echo on

wc=pi/3;

M=20;

MH=M/2;

n=(0:M);

for k=1:M+1

hlp(k)=(wc/pi)*clsinc(wc*(n(k)-MH));

end

pause;

(Fhlp,w) = clwfour(hlp,128);pause;

WRect = (ones(1:M+1));pause;

(FRect,w) = clwfour(WRect,128);pause;

WBart = 1 - abs((n-MH)./MH);pause;

(FBart,w) = clwfour(WBart,128);pause;

WHann = .5 - 0.5*cos(2*pi*n/M);pause;

(FHann,w) = clwfour(WHann,128);pause;

WHamm = 0.54 - .46*cos(2*pi*n/M);pause;

(FHamm,w) = clwfour(WHamm,128);pause;

WBlack = 0.42 - 0.5*cos(2*pi*n/M) + 0.08*cos(4*pi*n/M);pause;

(FBlack,w) = clwfour(WBlack,128);pause;

clg;

plot(n,WRect,n,WBart,n,WHann);pause;

plot(n,WHann,n,WHamm,n,WBlack);pause;

plot(w,abs(FRect),w,abs(FBart),w,abs(FHann));pause;

plot(w,abs(FHann),w,abs(FHamm),w,abs(FBlack));pause;

plot(w,20*log10(abs(FRect)),w,20*log10(abs(FBart))...

,w,20*log10(abs(FHann)));pause;

plot(w,20*log10(abs(FHann)),w,20*log10(abs(FHamm))...

,w,20*log10(abs(FBlack)));pause;

hHann=hlp.*WHann;pause;

(FhHann,w)=clwfour(hHann,128);pause;

clg;plot(w,abs(Fhlp),w,abs(FhHann));pause;

plot(w,20*log10(abs(Fhlp)),w,20*log10(abs(FhHann)));pause;

hBart=hlp.*WBart;pause;

(FhBart,w)=clwfour(hBart,128);pause;

clg;plot(w,abs(Fhlp),w,abs(FhBart));pause;

plot(w,20*log10(abs(Fhlp)),w,20*log10(abs(FhBart)));pause;

B=0;

Den=bessel(0,j*B);

for k=1:M+1;

hKais0(k)=bessel(0,j*B*(1-((n(k)-MH)/MH)^2)^.5)/Den;

end

(FKais0,w) = clwfour(hKais0,128);pause;

B=3;

Den=bessel(0,j*B);

for k=1:M+1;

hKais3(k)=bessel(0,j*B*(1-((n(k)-MH)/MH)^2)^.5)/Den;

end

(FKais3,w) = clwfour(hKais3,128);pause;

B=6;

Den=bessel(0,j*B);

for k=1:M+1;

hKais6(k)=bessel(0,j*B*(1-((n(k)-MH)/MH)^2)^.5)/Den;

end

(FKais6,w) = clwfour(hKais6,128);pause;

plot(n,hKais0,n,hKais3,n,hKais6);pause;

plot(w,abs(FKais0),w,abs(FKais3)...

,w,abs(FKais6));pause;

plot(w,20*log10(abs(FKais0)),w,20*log10(abs(FKais3))...

,w,20*log10(abs(FKais6)));pause;

T=1;

for k=1:MH

arg=pi*(n(k)-MH);

hdiff(k)=(arg*cos(arg)-sin(arg))/(pi*(n(k)-MH)*(n(k)-MH)*T);

hdiff(M+2-k)=-hdiff(k);

end;

hdiff(MH+1)=0.0;

(Fdiff,w) = clwfour(hdiff,128);pause;

hdiffhann=hdiff.*WHann;pause;

(FdiffHann,w)=clwfour(hdiffhann,128);pause;

FREQUENCY SAMPLING TECHNIQUES

From DFT work, we know that we can find M+1 values for h(n) for any M+1 values of H(k) that satisfy H(M+1-k) = H*(k) with a cosine series expansion of the form:

(Note: just like Chapter 7 notes with N instead of M+1).

For now, we won't worry about a possible "singleton" at a k value which might sample at w=p.

To simplify the 2 thing, we can rewrite the form equation as:

where

0<k<M/2

We would get the same values of h(n) if we just took the ifft, however the formula gives us the big advantage of making non-integer time shifts. We will demonstrate that a desirable form for an appropriate selection of LPF frequency samples is:

The half sample delay in the time domain (e ^jw/2 in the frequency domain) allows us to preserve the h(n) = h(M+1 - n) symmetry.

Example: M=16, pass to 3p/8, no alternating sign.

BLOW UP SHOWING THAT CRITERIA WERE MET IN THE FREQUENCY DOMAIN

CRITERIA STILL MET

ADD TRANSITION POINT TO DECREASE RIPPLING

BLOW-UP OF SPECS - of course the transition band has been increased.

If this is not satisfactory, we would need to increase N to make it sharper.

% clfirfs

% WHY WE NEED ALTERNATING SIGNS

wo=2*pi/16;

k=[0:15];

k64=[0:63];

n=k;

x1=exp(j*wo*2*n);

X1=fft(x1);samp(abs(X1));pause;

X1e=fft(x1,64);

plot(k64,real(X1e));pause;

x2=exp(j*wo*3*n);

X2e=fft(x2,64);

plot(k64,real(X2e));pause;

plot(k64,abs(X1e+X2e),k64,abs(X1e+X2e),'o');pause;

plot(k64,abs(X1e-X2e),k64,abs(X1e-X2e),'o');pause;

x3=exp(j*wo*4*n);

X3e=fft(x3,64);

plot(k64,real(X3e));pause;

plot(k64,abs(X1e+X2e+X3e));pause;

plot(k64,abs(X1e-X2e+X3e));pause;

% IDEAL LOW PASS, ZERO PHASES

Hk=[1 1 1 1 zeros(1,9) 1 1 1];

samp(k,Hk);pause;

wk=2*pi/16*k;

h=ifft(Hk);

samp(k,h);pause;

w=2*pi/64*k64;

Hw=clpfour(w,k,h);pause;

clg;plot(w,abs(Hw),wk,abs(Hk),'o');pause;

clg;plot(w,angle(Hw),wk,angle(Hk),'o');pause;

% IDEAL LOW PASS, ALTERNATING SIGNS

% LIKE DELAY OF M/2 IN THE TIME DOMAIN

Hka=[1 -1 1 -1 zeros(1,9) -1 1 -1];

ha=ifft(Hka);

clg;samp(k,ha);pause;

Hwa=clpfour(w,k,ha);pause;

clg;plot(w,abs(Hwa),wk,abs(Hka),'o');pause;

clg;plot(w,angle(Hwa),wk,angle(Hka),'o');pause;

%KEEP ALTERNATING SIGNS BUT ADD LINEAR PHASE SHIFT

wo=2*pi/16;

np2=n+.5;

haf=1/16*(1-2*cos(wo*np2)+2*cos(2*wo*np2)-2*cos(3*wo*np2));

clg;samp(n,haf);pause;

plot(n,ha,n,haf);pause;

Hwaf=clpfour(w,k,haf);pause;

clg;plot(w,abs(Hwaf),wk,abs(Hka));pause;

clg;plot(w,angle(Hwaf),wk,angle(Hka));pause;

% ADD TRANSITION POINT

Hkat=[1 -1 1 -1 .38 zeros(1,7) .38 -1 1 -1];

samp(abs(Hkat));pause;

hat=ifft(Hkat);

samp(k,hat);pause;

Hwat=clpfour(w,k,hat);pause;

clg;plot(w,abs(Hwat),wk,abs(Hkat),'o');pause;

clg;plot(w,angle(Hwat),wk,angle(Hkat),'o');pause;

The reduced value at P needed to effect a transition between the pass and stop bands, results in a 40 dB attenuation in the stop band and gives us a way to determine the resolution needed in the frequency domain:

This N may need to be increased if it doesn't match the passband specifications.

Book example: w_c = 0.08p, w_r = 0.20p, A > 40 dB

However, this is not sufficient because the stop band begins at .18p and the passband is only guaranteed up to .06p, i.e., the sampled points are at (in rads: [0 .06p .12p .18p .24p ...], so H = [1 -1 .38 0 0 ...] isn't good for the passband and H=[1 -1 1 -.38 0 ...] isn't good the stop band.

therefore a better resolution is:

EQUIRIPPLE DESIGNS PRINCIPLE

To eliminate dealing with phase problems, we go back to the
h(n) = h(-n) symmetry, i.e.

Then we use the trigonometric idea that:

cos nw can always be written as a polynomial in cos w

e.g., cos 2w	= cos2 w - sin 2w = cos2 w - (1- cos2 w)
	= -1 + 2 cos2 w

which has M/2 - 1 possible zeros for the polynomial part

+ two possible zeros at w = 0 and p

(the zeros give us the frequencies at which the ripples have max and min)

+ two points to be defined at w_c and w_r

giving us M/2 + 3 points to curve fit to polynomials.

Different iterative techniques have been used to do these fits.

The M value needed for a particular job is usually better then that used for other techniques. We looked at REMEZ in the first part of the book.

FOR THE EXAMPLE WE JUST DID:

N=15;

F = [0 3/8 4/8 1]; M = [1 1 0 0];

h=remez(N,F,M);

samp(h);pause;

H=fft(h,64);

plot([0:63],abs(H));pause;

Note: 3p/8 is at k=12 and 4p/8 is at k=16

OTHER TYPES OF FILTERS

What if we want something besides lowpass?

H'_HP(w) = 1 - H'_LP(w) -----> h _HP(n) = d(n-M/2) - h_LP(n)

H'_BP(w) = H'_LP2(w) - H'_LP1(w) -----> h_BP(n) = h_LP1(n) - h_LP2(n)

H'_BS(w) = 1 - H'_BP(w) -----> h_BS(n) = d(n-M/2) - h_LP2(n) + h_LP1(n)

Note: 1 is really exp(-j* w*M/2) since we need the delay for symmetry.

demo: clfirtr.m---patterned after problem 9.2

hlp1 = cllpdfir(pi/4,20);

hlp2 = cllpdfir(3*pi/4,20);

hhann=[0 hanning(39)' 0];

hlp1=hhann.*hlp1;

hlp2=hhann.*hlp2;

clwfour(hlp1,128);pause;

clwfour(hlp2,128);pause;

hhp=[zeros(1,20) 1 zeros(1,20)] - hlp1;

hbp=hlp2 - hlp1;

hbs=[zeros(1,20) 1 zeros(1,20)] - hbp;

clwfour(hhp,128);pause;

clwfour(hbp,128);pause;

clwfour(hbs,128);pause;

function [x] = cllpdfir(w,M)

for k=0:2*M;

x(k+1)=clsinc(w*(k-M))*w/pi; end;

LOWPASS hlp1(n) WITH w_d AT p/4

LOWPASS hlp2(n) WITH w_d AT 3p/4

d (n-20) - hlp1(n) =HIGHPASS WITH w_d AT p/4

hlp2(n)- hlp1(n) =hbp(n) = BANDPASS BETWEEN AT p/4 and 3p/4

d (n-20) - hlp2(n) + hlp1(n)= d (n-20) - hbp(n)=hbs(n)

=BANDSTOP BETWEEN AT p/4 and 3p/4