Transcript CS3291 : Digital Signal Processing '04-05

CS3291 : Digital Signal Processing '04-05
Section 3 : Discrete-time LTI systems
3.1
Introduction:
A discrete time system takes discrete time input signal { x[ n ] },
& produces an output signal { y[ n ] }.
input
{x[n]}
Oct'04
SYSTEM
CS3291
output
{y[n]}
1
• Implemented by general purpose computer, microcomputer or
dedicated digital hardware capable of carrying out arithmetic
operations on samples of {x[n]} and {y[n]}.
• {x[n]} is sequence whose value at t=nT is x[n].
• Similarly for {y[n]}.
• T is sampling interval in seconds.
• 1/T is sampling frequency in Hz.
• {x[n-N]} is sequence whose value at t=nT is x[n-N].
{x[n-N]}
is
{x[n]}
with
every
sample
delayed
by N sampling intervals.
• Many types of discrete time system can be considered, e.g.
Oct'04
CS3291
2
(i) Discrete time "amplifier" with output:
y[n] = A . x[n].
•
Described by "difference equation": y[n] = A x[n].
•
Represented in diagram form by a "signal flow graph” .
x[n]
Oct'04
A
y[n]
CS3291
3
(ii) Processing system whose output at t=nT is obtained by
weighting & summing present & previous input samples:
y[n] = A0 x[n] + A1 x[n-1] + A2 x[n-2] + A3 x[n-3] + A4x[n-4]
‘Non-recursive difference equatn’ with signal flow graph below.
Boxes marked " z -1 " produce a delay of one sampling interval.
x[n]
z-1
A0
z-1
A1
z-1
A2
z-1
A3
A4
y[n]
Oct'04
CS3291
4
(iii) System whose output y[n] at t = nT is calculated according to
the following recursive difference equation:
y[n] = A 0 x[n] - B 1 y[n-1]
whose signal flow graph is given below.
• Recursive means that previous values of y[n] as well as
present & previous values of x[n] are used to calculate y[n].
x[n]
A0
y[n]
B1
Oct'04
CS3291
z-1
5
(iv) A system whose output at t=nT is:
y[n] = (x[n]) 2
as represented below.
y[n]
x[n]
Oct'04
CS3291
6
3.2. LTI Systems:
To be LTI, a discrete-time system must satisfy:
(I) Linearity ( Superposition ):
•Given any two discrete time signals {x 1 [n]} & {x 2 [n]},
if {x 1 [n]}  {y 1 [n]} & {x 2 [n]}  {y 2 [n]}
then for any values of k 1 and k 2 ,
response to k 1{x 1[n]} + k 2{x 2[n]}
must be k 1{y1[n]} + k 2 {y 2 [n]} .
•To multiply a sequence by k, multiply each element by k,
k{x[n]} = {k x[n]}.
•To add two sequences together, add corresponding samples,
{x[n]} + {y[n]} = {x[n] + y[n]}.)
Oct'04
CS3291
7
(II) Time-invariance :
• Given any discrete time signal {x[n]},
if response to {x[n]} is {y[n]},
response to {x[n-N]} must be {y[n-N]} for any N.
• Delaying input by N samples only delays output by N samples.
• Examples (i), (ii) and (iii) are LTI whereas (iv) is not LTI .
Oct'04
CS3291
8
3.3. Impulse-response:
Useful to consider response of LTI systems to a discrete time unit
impulse , or in short an impulse denoted by {d [n] } with:
1 : n  0
d[n]  
0 : n  0
d[n]
d[n-4]
1
1
...
-3 -2 -1 0 1
Oct'04
2
3
4
5
...
n
n
...
CS3291
-3 -2 -1 0 1
2
3
4
5
...
9
• {d[n-N]} is delayed impulse where the only non-zero sample
occurs at n=N rather than at n=0.
• When input is {d[n]}, output is impulse response {h[n]}.
• If impulse-response of an LTI system is known, its response
to any other input signal may be predicted.
Oct'04
CS3291
10
3.4. Computing an impulse response:
• Consider a discrete time system with signal-flow graph below
with A1, A2, A3, A4, A5 set to specific constants.
• Notice the labels X1, X2, etc.
• Realised by microcomputer running a program with flow-diagram
& high-level language implementation as listed next.
X1
z-1
x[n]
A1
X2
z-1
A2
X3
z-1
A3
X4
z-1
A4
X5
A5
Y
y[n]
Oct'04
CS3291
11
clear all;
A1=1; A2=2; A3=3; A4=-4; A5=5;
Set A1,A2,….
x1=0; x2=0; x3=0; x4=0; x5=0;
Zero X1, X2, X3, X4,X5
while ~isempty(X)
X = input( 'X = '); x1 = X;
INPUT X
Y= A1*x1 + A2*x2 + A3*x3
Y = A1*X + A2*X2 + A3*X3 +...
+ A4*x4 + A5*x5 ;
disp(['
OUTPUT Y
X5 = X4; X4 = X3; X3 = X2; X2 =X1
Y = ' num2str(Y)]);
x5 = x4 ; x4 = x3;
x3 = x2 ; x2 =x1;
end;
% MATLAB program
Oct'04
CS3291
12
Here is another version using MATLAB arrays:
A = [1 2 3 -4 5 ] ;
x = [0 0 0 0 0 ] ;
while 1
X = input( 'X = '); x(1)=X;
Y=A(1) * X;
for k = 2 : 5
Y = Y + A(k)*x(k);
end;
disp(['
Y = ' num2str(Y)]);
for k=5:-1:2
x(k) = x(k-1);
end;
end;
Oct'04
CS3291
13
And an even more efficient version
A = [1 2 3 -4 5 ] ;
x = [0 0 0 0 0 ] ;
while 1
x(1) = input( 'X = ');
Y = A(1)*x(1);
for k = 5 : -1: 2
Y = Y + A(k)*x(k);
x(k) = x(k-1);
end;
disp(['
Y = ' num2str(Y)]);
end;
Oct'04
CS3291
14
• These can be run in MATLAB
• The ‘while 1’ statement initiates an infinite loop.
• Program will run for ever
or until interrupted by ‘CONTROL+C’
Oct'04
CS3291
15
• Impulse-response for this example may be deduced by running
program & entering values for X on the keyboard :
0, 0, 0, 1, 0, 0, 0, 0,....
• Sequence of output samples printed out will be :
0, 0, 0, A1, A2, A3, A4, A5, 0, 0, ....
• Impulse-response can also be obtained by tabulation (later).
• Output must be zero until input becomes equal to 1 at n=0
• Impulse response is:
{..., 0, ..., 0, A1, A2, A3, A4, A5, 0, 0, ... ,0,...}
where the sample at n=0 is underlined.
• Only five non-zero output samples are observed.
• This is a " finite impulse-response ” (FIR).
Oct'04
CS3291
16
Exercise 3.1:
Calculate impulse-responses for difference equations:
(i) y[n] = x[n] + 2 x[n-1] + 3 x[n-2] - 4 x[n-3] + 5 x[n-4]
(ii) y[n] = 4 x[n] - 0.5 y[n-1]
Solutions:
• If computer not available we may use a tabulation.
• Difference eqn (i) will produce a finite impulse-response.
• Difference eqn (ii) produces infinite response whose samples
gradually reduce in amplitude but never quite become zero.
Oct'04
CS3291
17
(i) y[n] = x[n] + 2 x[n-1] + 3 x[n-2] - 4 x[n-3] + 5 x[n-4]
n
-1
0
1
2
3
4
5
:
x[n]
0
1
0
0
0
0
0
:
x[n-1]
0
0
1
0
0
0
0
:
x[n-2]
0
0
0
1
0
0
0
:
x[n-3]
0
0
0
0
1
0
0
:
x[n-4]
0
0
0
0
0
1
0
:
y[n]
0
1
2
3
-4
5
0
:
Solution is:- {.. 0, .., 0, 1, 2, 3, -4, 5, 0, .., 0, ...}
Oct'04
CS3291
18
Difference equation (ii): y[n] = 4 x[n] - 0.5 y[n-1]
n
0
1
2
3
4
5
:
x[n]
1
0
0
0
0
0
:
y[n-1]
0
4
-2
1
-0.5
0.25
:
y[n]
4
-2
1
-0.5
0.25
-0.125
:
Solution is: {.., 0, .., 0, 4, -2, 1, -0.5, 0.25, -0.125, ...)
Oct'04
CS3291
19
3.5. Digital filters:
•A "digital filter” is a digitally implemented LTI discrete time system
governed by a difference equation of finite order; e.g. :
(i) y[n] = x[n] + 2 x[n-1] + 3 x[n-2] - 4 x[n-3] + 5 x[n-4]
(ii) y[n] = 4 x[n] - 0.5 y[n-1]
• Difference equation (i) is "non-recursive"
& produces a finite impulse response (FIR).
• Difference equation (ii) is " recursive " .
• Impulse-response of a recursive difference equation can have an
infinite number of non-zero terms.
In this case it is an infinite impulse-response (IIR).
Oct'04
CS3291
20
• Next slide has yet another way of programming difference
equations (i.e. digital filters) in MATLAB.
• But this time, we use a MATLAB function and operate on a
whole array of music at once.
• You can find the music file ‘cap4th.wav’ in
•www.cs.man.ac.uk/~barry/mydocs/CS3291/MATLAB
Oct'04
CS3291
21
MATLAB function ‘filter’ used for non-recursive diffce eqn
y = filter(A, 1, x) filters signal in array x to create array y.
For each sample x(n) of array x:
y(n) = A(1)*x(n) + A(2)*x(n-1) + ... + A(N+1)*x(n-N)
Example (i):
y[n] = x[n] + 2 x[n-1] + 3 x[n-2] - 4 x[n-3] + 5 x[n-4]
[x, fs, nbits] = wavread ('cap4th.wav');
sound(x,fs,nbits);
pause;
y = filter ( [1 2 3 -4 5], 1, x);
sound(y,fs,nbits);
Oct'04
CS3291
22
Use of ‘filter’ for recursive difference equatn
y = filter(A, B, x) filters signal in array x to create array y
using a recursive difference equation.
Example (ii): y[n] = 4 x[n] - 0.5 y[n-1]
[x, fs, nbits] = wavread ('cap4th.wav');
sound(x,fs,nbits);
pause;
y = filter ( [4], [1 0.5], x);
sound(y,fs,nbits);
Oct'04
CS3291
23
3.6. Discrete time convolution:
If impulse-response of an LTI system is {h[n]} its response to any
input {x[n]} is an output {y[n] } whose samples are given by the
following "convolution" formulae:
y[n] 

 h[m] x[n  m]
or
y[n] 

 x[k] h[n  k]
k  
m  
Consider the response of a system with impulse response:
{h[n]} = { ..., 0,..., 0, 1, 2, 3, -4, 5, 0, .....0, .... }
to {x[n]} = { ... 0, ... , 0 , 1, 2, 3, 0, ..., 0, ....}
y[n] 
4
 h[m] x[n  m]
 x[n]  2x[n  1]  3x[n  2]  4x[n  3]  5x[n  4]
m 0
Oct'04
CS3291
24
• Result is not too surprising.
• It is the difference equation for an LTI system whose impulse
response is {.., 0, .., 0, 1, 2, 3, -4, 5, 0, .., 0, ..}.
• Program discussed earlier implements this difference equation,
•We could complete the example by running it, entering samples of
{ x[n] }, and noting the output sequence.
• Alternatively, we could use tabulation as follows:
Oct'04
CS3291
25
y[n] = x[n] + 2x[n-1] + 3x[n-2] - 4x[n-3] + 5x[n-4]
n
x[n] x[n-1] x[n-2] x[n-3] x[n-4] y[n]
:
:
:
:
:
:
:
-1
0
0
0
0
0
0
0
1
0
0
0
0
1
1
2
1
0
0
0
4
2
3
2
1
0
0
10
3
0
3
2
1
0
8
4
0
0
3
2
1
6
5
0
0
0
3
2
-2
6
0
0
0
0
3
15
7
0
0
0
0
0
0
:
:
:
:
:
:
:
{y[n]} = { .... 0, ....., 0, 1, 4, 10, 8, 6, -2, 15, 0, ...., 0, ....}
Oct'04
CS3291
26
3.7. Proof of discrete time convolution formulae:
Only use defns of linearity, time-invariance & impulse-response.
Since d[n-m] is non-zero only at n = m, given any sequence {x[n]},
x[n] 

 x[m]d[n  m]
m  
 {x[n]} 

 x[m]{d[n  m]}
m  
{x[n]} is sum of infinite number of delayed impulses{d[n-m]}
each multiplied by a single element, x[m].
Response to {d[n-m]} is {h[n-m]} for any value of m.
 response to {x[n]} is :
{y[n]} 

 x[m]{h[n  m]}
i.e. y[n] 
m  
Oct'04

 x[m]h[n  m]
for all n
m  
CS3291
27
• Replacing n-m by k gives the alternative formula.
• Study the graphical explanation of this proof in Section 3.17.
x[n]
{x[n]} = { ..., 0, ..., 0, 1, 2, 3, 0, 0, ..., 0, ...}
3
2
n
h[n]
{h[n]} = { ..., 0, ..., 0, 1, 2, 1, 2, 0, ..., 0...}
2
n
Oct'04
CS3291
28
3.17 Discrete time convolution graphically
x[n]
{h[n]} = { ..., 0, ..., 0, 1, 2, 1, 2, 0, ..., 0,...}
3
2
{x[n]} = { ..., 0, ..., 0, 1, 2, 3, 0, 0, ..., 0, ...}
h[n]
n
d[n]
2
n
n
2d[n-1]
2h[n-1]
4
2
2
n
n
3d[n-2]
3h[n-2]
6
3
3
n
Oct'04
CS3291
n
29
y[n]
10
8
6
4
2
n
Express:
{x[n]} = {d[n]} + 2 {d[n-1]} + 3 {d[n-2]}
Response is:
{y[n]} = {h[n]} + 2 {h[n-1} + 3 {h[n-2]}
Oct'04
CS3291
30
x[n]
{h[n]} = { ..., 0, ..., 0, 1, 2, 1, 2, 0, ..., 0,...}
{x[n]} = { ..., 0, ..., 0, 1, 2, 3, 0, 0, ..., 0, ...}
3
2
n
{ d[n] }
{h[n]}
+ 2{ d[n-1 ]}
+ 2 { h[n-1]}
+ 3 { d[n-2] }
+ 3 { h[n-2] }
y[n]
10
8
6
4
2
Oct'04
CS3291
n
31
d[n]
h[n]
2
n
n
2d[n-1]
2h[n-1]
4
2
2
n
n
3d[n-2]
3h[n-2]
6
3
3
n
Oct'04
n
CS3291
32
3.8. Stability:
An LTI system is stable if its impulse-response {h[n]} satisfies:

 h[n]
is finite
n  
This means that {h[n]} must be either an FIR or an IIR whose
samples decay towards zero as n  .
3.9. Causality:
An LTI system operating in real time must be "causal " which
means that its impulse-response {h[n]} must satisfy:
h[n] = 0 for n < 0.
Non-causal system would need “crystal ball ” to predict future.
Oct'04
CS3291
33
3.10. Relative Frequency
• Study effect of digital filters on sine-waves of different frequencies.
• Discrete time sine-wave obtained by sampling A cos(t + ).
• If sampling period is T seconds, resulting discrete time signal is :
{A cos(Tn   )}  {A cos( n   )}
•  = T is “ relative" frequency of sampled sine-wave.
• Units of  are 'radians per sample'.
• To convert  back to true frequency (radians/s ) multiply by 1/T.
• Note: 1/T = sampling frequency, f S, in samples/second (Hz).
‘radians / sample’ times ‘samples / second’ = ‘radians / second’
•Normally process analogue signals in range 0 to f S /2 ( = 1/(2T) ),
• Restricts  to the range 0 to .
Oct'04
CS3291
34
cos( t )
-4T
-3T
=2 / (8T)
3T
-T
t
4T
T
cos( T n )
 = 2 / 8
-4
-3
3
-1
Oct'04
n
4
1
CS3291
35
Table gives values of  & corresponding true frequencies :Relative frequency 
True frequency
(radians/sample)
(radians/s) (Hz)
0
0
0
/6
 fS/6
fS/12
/4
 fS/4
fs/8
/3
 fS/3
fs/6
/2
 fS/2
fs/4
2/3
2fS/3
fs/3

fS
fS/2
‘radians / sample’  ‘samples / second’ = ‘radians / second’
Oct'04
CS3291
36
3.11. Relative frequency response
Response to sampled sine-waves calculated by defining:
{e jn }  {cos( n)  j sin( n)}
If this is applied to a system with impulse response {h[n]}, output
would be {y[n]} with :
y[n] 


h[m]e
j ( n  m )

m  
e
jn

jn  jm
h
[
m
]
e
e

m  

jn
h
[
m
]
e

 H (e j )e jn
m  
j
where H (e ) 

 h[m]e
m  
Oct'04
CS3291
 jm


 jn
h
[
n
]
e
(DTFT)

n  
37
cos( (2/8) n )
-4
-3
3
-1
n
4
1
sin( (2/8) n)
-2
n
-1
1
Oct'04
CS3291
3
38
• H( e j  ) is the " relative frequency response " of LTI system.
• Complex number for any value of .
• It is DTFT of { h[n] }.
• Note similarity to the analogue Fourier transform.
j
H (e ) 

 h[n]e
 jn
(DTFT)
n  

H a ( j )   ha (t )e
 jt

Oct'04
CS3291
(FT)
39
• If input is {e
j  n
}, output is
same sequence with each
element multiplied by H( e j  ).
{e j  n }
{ e j  n H(e j  n ) }
LTI
{h[n]}
j
H (e ) 

 jn
h
[
n
]
e

n  
•When { h[n] } is real, H( e -j  ) = H*( e j  ).
•Exercise: Prove this.
Oct'04
CS3291
40
• Remember the analogue case:-
e j t
e j t Ha(j )
LTI
ha(t)

H a ( j )   ha (t )e jt

• When h(t) is real, Ha(j ) = Ha*(j )
Oct'04
CS3291
41
3.12. Gain and phase responses
• Expressing H( e j  ) = G() e j ()
• G() = |H( e j  )| is "gain ”
• () = arg ( H( e j  ) ) is "phase lead”.
• Both vary with .
• If input is {A cos(n)},
•
output is: { G()A cos(n + ()) }
Exercise: Prove this
•When input is sampled sine-wave of rel frequency ,
output is sine-wave of same frequency ,
but with amplitude scaled by G() & phase increased by ().
Oct'04
CS3291
42
Gain & phase response graphs again
-() 20log10[G()] dB
-()
G() in dB
0
Oct'04
/4
/2
3/4
CS3291


43
•Common to plot graphs of G() & () against .
•G() often converted to dBs by calculating 20 log10( G() ).
•Restrict  to lie in range 0 to  & adopt a linear horizontal scale.
•Example 3.2:
Calculate gain & phase response of FIR digital filter in Figure 3.3
where A1 = 1, A2 = 2, A3 = 3, A4 = -4, A5 = 5 .
•Solution:
Impulse response is: {.., 0, .., 0, 1, 2, 3, -4, 5 ,0, .., 0,..}.
By the formula established above,
H( e j  ) = 1 + 2 e - j  + 3 e -2 j  - 4 e -3 j  + 5 e - 4 j 
To obtain the gain and phase at any  take modulus and phase.
Best done by computer, either by writing & running a simple
program, or by running a standard package.
A simple program is given at the end of these notes.
Oct'04
CS3291
44
• Can use a spread-sheet such as excel (see lecture notes)
• Another way of calculating gain & phase responses for
H( e j  ) = 1 + 2 e - j  + 3 e -2 j  - 4 e -3 j  + 5 e - 4 j 
is to start MATLAB and type:
A = [ 1 2 3 -4 5 ];
freqz( A );
• See later section 3.18 for graphs produced.
Oct'04
CS3291
45
• Produces graphs of gain & phase responses of H( ej ).
• Frequency scale normalised to fS/2 & labelled 0 to 1
instead of 0 to .
• See Section 3.18
• In next section we calculate coeffs of an FIR filter to achieve a
particular type of gain & phase response.
• To calculate gain & phase responses for an IIR filter would be
difficult by method used above because impulse-response has
infinite number of non-zero terms.
•There is an easier method in a later section.
Oct'04
CS3291
46
Phase delay
• Assume input is {A cos(n)},
• Output is: { G() A cos (  n
+
() ) }
= { G() A cos (  [ n + () /  ] ) }
= { G() A cos (  [ n - (-() / ) ] ) }
•Phase lead () delays sine-wave by -()/ sampling intervals
• This is ‘phase-delay’.
Oct'04
CS3291
47
3.13. Linear phase response
• If -()/ remains constant for all values of ,
i.e. if -() = k for constant k,
system has ‘linear phase’ response.
• Graph of -() against  on linear scale would then be
straight line with slope k where k is "phase delay"
i.e. the delay measured in sampling intervals.
• This need not be an integer.
Oct'04
CS3291
48
Gain, phase & ‘phase-delay’ response graphs
-()/
-()
-()/
G() in dB
0
/4
/2
3/4


Not ‘linear-phase’.
Oct'04
CS3291
49
Gain, phase & ‘phase-delay’ response graphs again
-()/
-()
-()/
k
G() in dB
0
/4
/2
3/4


This is ‘linear-phase’ as -()/ is constant.
Oct'04
CS3291
50
• Linear phase systems of special interest since all sine-waves are
delayed by same number of sampling intervals.
• Input signal expressed as Fourier series will have all its
sinusoidal components delayed by same amount of time.
• In absence of amplitude changes, output signal will be exact
replica of input.
•LTI systems are not necessarily linear phase.
• A certain class of digital filter can be made linear phase.
Oct'04
CS3291
51
3.14. Inverse DTFT
• Frequency -response H( e j  ) is DTFT of {h[n]}.
• Inverse DTFT formula allows {h[n]} to be deduced from H(e j ) :
1
h[n] 
2



H (e j )e jn d
: - < n < 
• Observe similarities with analogue Fourier transform:
1
ha (t ) 
2



H a ( j )e jt d
: - < t < 
(i) (1/2) factor,
(ii) sign of jn (n replaces t).
(iii) variable of integration (d)
Oct'04
CS3291
52
Also note that with DTFT:
• Negative frequencies are again included.
• Range of integration is now - to  rather than - to .
• Forward transform is a summation
& inverse transform is an integral.
• This is because {h[n]} is a sequence
whereas H(e j  ) is function of the continuous variable .
Oct'04
CS3291
53
3.15 Problems:
1. Show that y[n]=(x[n]) 2 is not LTI.
2. If {h[n]}= { .. 0, .. , 0, 1, -1, 0, .. 0, .. }, calculate response to
{x[n]} = { .... 0, .., 0, 1, 2, 3, 4, 0, .., 0, ... }.
3.. Produce a signal flow graph for each of the following difference equations:
(i) y[n] = x[n] + x[n-1] (ii) y[n] = x[n] - y[n-1]
4. For each difference equation in question 3, determine the impulse response,
& deduce whether the system it represents is stable and/or causal.
.5. Calculate, by tabulation, the output from difference equation (ii) in question 3
when the input is the impulse response of difference equation (i).
6. If fS =8000 Hz, what true frequency corresponds to /5 radians/sample?
7. Sketch the gain & phase responses of the system referred to in question 2.
(You don't need a computer program.) Is it a linear phase system?
8. Calculate the impulse response for y[n] = 4x[n] - 2y[n-1]. Is it stable?
9. Show that input {cos(n)} produces an output {G()cos(n+())}.
10. For y[n]=x[n]+2x[n-1]+3x[n-2]+2x[n-3]+x[n-4], sketch phase response and
comment on its significance. Show that () = -2.
Oct'04
CS3291
54
3.16. Progam for gain & phase responses of FIR digital filter.
(using complex arithmetic)
fs = input('enter fs:'); M = input('Enter order:');
disp('Enter coefficients:');
for N = 0:M C(N+1)=input('coeff = '); end;
disp('
FREQ(Hz)
GAIN(dB.)
PHASE LEAD(deg)');
K = 0; F =0.0; % F is frequency in Hz}
while F < fs/2
H=0;
W =2*pi*F/fs;
for N=0:M
H=H+C(N+1)*exp(-j*W*N); end;
MD=abs(H); PH = angle(H); GAIN(K+1)=-999;
if MD>0 GAIN(K+1)=20*log10(MD);end;
PHAS(K+1)=PH*180.0/pi;
disp([num2str(F) ' 'num2str(GAIN(K+1)) '
F=F+fs/200; K=K+1;
'num2str(PHAS(K+1))]);
end; %while loop
Oct'04
CS3291
55
3.16. Progam for gain & phase responses of FIR digital filter.
(without complex arithmetic)
fs = input('enter fs:'); M = input('Enter order:');
disp('Enter coefficients:');
for N = 0:M H(N+1)=input('coeff = '); end;
disp(' FREQ(Hz) GAIN(dB.) PHASE LEAD(deg)');
K = 0; F =0.0; % F is frequency in Hz}
while F < fs/2
Hre =0.0; Him =0.0; W =2*pi*F/fs;
for N=0:M
Hre=Hre+H(N+1)*cos(W*N); Him=Him-H(N+1)*sin(W*N); end;
MD=sqrt(Hre*Hre+Him*Him);
PH = atan2(Him,Hre); GAIN(K+1)=-999;
if MD>0 GAIN(K+1)=20*log10(MD);end;
PHAS(K+1)=PH*180.0/pi;
disp([num2str(F) ' 'num2str(GAIN(K+1)) ' 'num2str(PHAS(K+1))]);
F=F+fs/200; K=K+1;
end; %while loop
Oct'04
CS3291
56
Sample run:FIR FREQUENCY RESPONSE
Enter sampling frequency(Hz): 1000
Enter order of FIR filter: 3
Enter coefficients:H[0] = 1 H[1] = 2 H[2] = 1 H[3] = 2
Output obtained:Freq(Hz) Gain(DB)
0
15.56
5
15.56
10
15.56
:
:
:
500
6.01
Oct'04
Phase(deg)
0.0
-3.0
-6.0
-180.0 }
CS3291
57
3.18. Frequency response of FIR digital filter by MATLAB
• Given discrete time LTI system with difference equation:
y[n] = x[n] + 2 x[n-1] + 3 x[n-2] - 4 x[n-3] + 5 x[n-4]
• Its impulse response is: { ... , 0, ... , 0, 1 , 2, 3, -4, 5, 0, ..., 0, ... }
• Its frequency response is:
H( e j  ) = 1 + 2 e - j  + 3 e -2 j  - 4 e -3 j  + 5 e - 4 j 
• To calculate its gain and phase responses, type:
freqz ( [ 1 2 3 -4 5 ]) ;
• Resulting graphs have a normalised frequency scale
with “1” corresponding to  =  radians/sample;
i.e. one half of the sampling frequency.
Oct'04
CS3291
58
Oct'04
CS3291
59
Oct'04
CS3291
60