Transcript Quiz3_Review_FilterDesign

Correlation and Spectral
Analysis
Application 4
Review of covariance
For independen
t random variables , x and y .
z  x-y forms a new random variable
z with vari
 z   x   y  cov( x , y )
2
2
2
cov( x , y )  ( x  x )( y  y  0 (independe
nce)
ance :
Autocorrelation (Autocovariance)
Recall that
2
σn 
1
N
N

(nq )
2
q0
for noise process n .
In general,
Rn ( j) 
1
N
N
n
q0
q
n q  j , where j  0 captures covariance
Noise Power
Rn ( j) 
Sk 
N
1
n
N
1
N
nq j
q0
2
N
n
j
e
 2i
jk
N
 DFT [ R n ( j )]
j0
For white
noise,
Sk 
2
Nk
q

2
by the DC value theorem
Zero-Mean Gaussian Noise
Power Spectrum
E{Pn(k)}  2 = 1.12 = Rn(0)
Auto-correlation
Rn ( j) 
1
N
N
n
q
nq j ,
q 1
>> for j = 1:256,
R(j) = sum(n.*circshift(n',j-1)');
end
Rn(0)  2 = 1.12
Window Selection: Hamming
y = filter(Hamming,1,n);
Hamming Filtered Power
Spectrum
White Noise Auto-Covariance
vs. Hamming Filtered Noise
Filtered
Noiseimage = imnoise(I,’gaussian’,0,10);
N_autocov = xcorr2(Noiseimage);
figure;imagesc(N_autocov/(128*128));colormap(gray);axis('image')
Image Noise Field
Autocovariance
Unfiltered
figure;imagesc(fftshift(abs(fft2(N_autocov/(128*128)))));colormap(gray);axis('image')
Image Noise Field
Power Spectrum
Filtered (wc = 0.6; order 20; Hamming Window)
N_autocov = xcorr2(Noiseimage_filtered);
figure;imagesc(N_autocov/(128*128));colormap(gray);axis('image')
Image Noise Field
Autocovariance
Filtered (wc = 0.6; order 20; Hamming Window)
N_autocov = xcorr2(Noiseimage_filtered);
figure;imagesc(N_autocov/(128*128));colormap(gray);axis('image')
Image Noise Field
Power Spectrum
fMRI Simulation
Windowing vs. Filtering
• “Window” applied in temporal or spatial domain
to reduce spectral leakage and ringing artifact
– Windows fall into a specialized set of functions
generally used for spectral analysis
• “Filter” applied to reduce noise, i.e. noise
matching, or to degrade or improve spatial
resolution
– Some cross-over: one method of filter design is the
“window” method which uses window functions for
frequency space modulating functions.
Windowing vs. Filtering
• Mathematically,
g ( x )  f ( x )  p ( x ) " Filter"
g ( t )  f ( t )  w ( t ) " Window"


G ( )  F ( )  P ( )
G ( f )  F ( f ) W ( f )
Filtering
MP 574
Outline
• Review of FIR/IIR Filters
– Z-transform
– Difference Equation
– Filter Design by Windowing
• Power Spectra
– Correlation and Convolution
– Example from Prof. Holden’s Notes
• Windowing and Spectral Estimation
• Weiner/Adaptive Filters
• Deconvolution
z-Transform as an Analysis Tool
• Sampled version (discrete version) of
the Laplace transform:
• z esT, where T is the sampling period.
• DFT and z-transform are related:
z = eiwT where s  eiwT

H (z) 

n  
h(n) z
n
ze
iw T
Laplace to z-Transform

b
F (s) 
e
 st
H (z) 
f ( t ) dt ,

h(n) z
n
n  
a
s    iw
Im(z)
Non-causal
signals
iw

unit circle
Re(z)
0
ws
Discrete FT
Continuous FT
z-Transform and Linear
Systems
• Stated more generally:
T{f(n)}
f(n)
h(n)

y (n) 
g(n)
 h(k ) x(n  k )
n0

H (z) 

h(n) z
n
, The inverse
z - transform
n0
H (z) 
N (z)
D (z)

ao z
N
 a1 z
N 1
 a2 z
N 2
 ...  a N
bo z
N
 b1 z
N 1
 b2 z
N 2
 ...  b N
Difference Equation
Implementation
• Shift theorem of z-transform:
bk x ( n  k )  bk z
k
X (z)
N
H (z) 
Y (z)
X (z)


bk z
k
k 0
M
1

ak z
k
k 0
N
y (n) 
b
k 0
M
k
x(n  k ) 
a
k 1
k
y (n  k )
Difference Equation
Implementation
• Shift theorem of z-transform:
bk x ( n  k )  bk z
k
X (z)
N
H (z) 
Y (z)
X (z)


bk z
k
k 0
M
1

ak z
k
k 0
N
y (n) 
b
k 0
M
k
x(n  k ) 
(n  k )
 a yFIR
k
k 1
FIR Coefficients and Impulse
Response
• FIR filter:
N
y (n) 
b
k
x(n  k )
k 0
b k  h ( k ), for FIR filters
FIR vs. IIR filters
• Finite impulse response (FIR) implies a
linear system that is always stable
– There are no poles
• Infinite impulse response (IIR) is only
stable if poles are inside the unit circle or
pole coincides with a zero.
IIR System
Im(z)
Zeros (o) at: -1, 2
Poles (x) at: 0.5±0.5j, 0.75
unit circle
x
o
Re(z)
x
x
o
IIR Stability
H1(z) 
H 2 (z) 
H 3(z) 
( z  1)( z  2 )
( z  0 . 5 j )( z  0 . 5 j )( z  0 . 75 )
( z  1)( z  2 )
( z  0 . 5 j )( z  0 . 5 j )( z  1 . 0 )
( z  1)( z  2 )
( z  0 . 5 j )( z  0 . 5 j )( z  1 . 1)
fvtool(B,A)
B = [1 -1 -2];
A=[1 -1.75 1.25 -0.375]
fvtool(B,A)
fvtool(B,A)
Unstable
B = [1 -1 -2];
A=[1 -1.75 1.25 -0.6]
Unstable
B = [1 -1 -2];
A=[1 -1.75 1.25 -0.6]
Finite impulse response (FIR)
B = [1 -1 -2];
A=[1]
Definition of Stability


k 0
h(k )  
FIR filter Design by Windowing
• Simply truncate IIR filter
• Rectangular Window:
 h d ( n ),
h(n)  
 0,
N1  n  N 2
Otherwise
h ( n )  hd ( n )  w ( n )
H (e
jw
)  H d (e
jw
) W (
Matlab: fdatool
filter() in Matlab
FILTER One-dimensional digital filter.
Y = FILTER(B,A,X) filters the data in vector X with the
filter described by vectors A and B to create the filtered
data Y. The filter is a "Direct Form II Transposed"
implementation of the standard difference equation:
a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
- a(2)*y(n-1) - ... - a(na+1)*y(n-na)
Exporting Filter Coefficients
Extension to 2D
• Radial Transform
– H(k)-> (H(k1)2+H(k2)2)1/2=T(k1,k2)
– See Matlab script on filter design using radial
transformation to 2D: Filter Design
• http://zoot.radiology.wisc.edu/~fains/Code/MP574_FilterDesign.m
• Parks-McClellan Transformation
– Step 1: Translate specifications of H(w1,w2) to H(w)
– Step 2: Design 1D filter H(w)
– Step 3: Map to 2D frequency space
cosw = - ½ + ½ cosw1 + ½ cosw2 + ½ cosw1 cosw2 = T(w1,w2)
- Step 4: determine h(n1,n2) by 2D FT.
Hamming Window Example
 2 
w ( n )  0 . 54  0 . 46 cos 
n ,
 N 
N /2
has form : H( w ) 

b ( n )(cos w )
n
N
2
(n N / 2)
, ... ,  1, 0 ,1, ... ,
N
2
 Type I, zero phase FIR filter.
n N / 2
Parks - McClellan
H( w 1 , w 2 )  H( w )
cos w  T ( w 1 , w 2 )
N /2


0 . 54  0 . 46 ( - ½  ½ cos w 1  ½ cos w 2  ½ cos w 1 cos w 2 )
n N / 2
n  N/2
Hamming Window Example
>> w1 = -pi:0.01:pi;
>> w2 = -pi:0.01:pi;
>> [W1,W2] = meshgrid(w1,w2);
>> H_2d = 0.54+0.46.*(-0.5+0.5.*cos(W1)+0.5.*cos(W2)+0.5.*cos(W1).*cos(W2));
>>figure;mesh(H_2d)
filter2()
2D FIR Filter Design,
Parks-McClellan
“firdemo”