Digital Filters

1
Digital Filters

• Part B Analysis, applications, and separability

2
Plotting filter frequency response(in
matlab/octave)
3
Impulse response for FIR filter

• Perform the z-transform (the discrete version of
  the Laplace transform) of h resulting H.
• Plot H on the unit circle. The magnitude of H
  (abs(H) or H) is amplitude and the angle of H
  (arg(H)) is the phase.
• Say we have a 3 point box filter
• h-1 h0 h1 1/3.

4
Matlab

• Matlab
• Commercial product
• The language of technical computing.
• MATLAB is a high-level language and interactive
  environment for numerical computation,
  visualization, and programming. Using MATLAB, you
  can analyze data, develop algorithms, and create
  models and applications. The language, tools,
  and built-in math functions enable you to explore
  multiple approaches and reach a solution faster
  than with spreadsheets or traditional programming
  languages, such as C/C or Java.
• You can use MATLAB for a range of applications,
  including signal processing and communications,
  image and video processing, control systems, test
  and measurement, computational finance, and
  computational biology. More than a million
  engineers and scientists in industry and academia
  use MATLAB, the language of technical computing.

5
GNU Octave

• GNU Octave is a high-level interpreted language,
  primarily intended for numerical computations.
  It provides capabilities for the numerical
  solution of linear and nonlinear problems, and
  for performing other numerical experiments. It
  also provides extensive graphics capabilities for
  data visualization and manipulation. Octave is
  normally used through its interactive command
  line interface, but it can also be used to write
  non-interactive programs. The Octave language is
  quite similar to Matlab so that most programs are
  easily portable.

6
(No Transcript)
7
amplitude

• w 0pi/100pi
• H 0.33exp(-1jw) 0.33exp(0jw)
  0.33exp(1jw)
• plot( w, abs(H) )
• grid
• xlabel( 'frequency \omega' )
• ylabel( 'amplitude' )

8
phase

• w 0pi/100pi
• H 0.33exp(-1jw) 0.33exp(0jw)
  0.33exp(1jw)
• plot( w, angle(H) )
• grid
• xlabel( 'frequency \omega' )
• ylabel( 'phase' )

9
amplitude (w/ more weight to center)

• w 0pi/100pi
• H 0.25exp(-1jw) 0.5exp(0jw)
  0.25exp(1jw)
• plot( w, abs(H) )
• grid
• xlabel( 'frequency \omega' )
• ylabel( 'amplitude' )

10
amplitude (highpass)

• w 0pi/100pi
• H -1exp(-1jw) 2exp(0jw)
  -1exp(1jw)
• plot( w, abs(H) )
• grid
• xlabel( 'frequency \omega' )
• ylabel( 'amplitude' )

11
IIR filter example
12
If we mistakenly do an in-place filter for the
homework

• yn 0.33 yn-1 0.33 xn 0.33 xn1
• w 0pi/100pi
• H (0.33exp(0jw) 0.33exp(1jw)) ./ (1 -
  0.33exp(-1jw))
• plot( w, abs(H) )
• grid
• xlabel( 'frequency \omega' )
• ylabel( 'amplitude' )

Note ./ is element-wise (pair-wise) division.
13
Application of filters

• Edge detection

14
Recall the first derivative test from calculus

• Let the function f be continuous on some interval
  (c-?,c?) containing the critical point c.
• 1. If f(x)gt0 for x in (c-?,c) and f(x)lt0 for x
  in (c,c?), then f has a local maximum at c.

15
Recall the first derivative test from calculus

• Let the function f be continuous on some interval
  (c-?,c?) containing the critical point c.
• 2. If f(x)lt0 for x in (c-?,c) and f(x)gt0 for x
  in (c,c?), then f has a local minimum at c.

16
Can we calculate the first derivative of an image?

• Let h1-1, 1.
• Then f(x) g(x) f(x) ? h1(x)
• where ? is cross correlation (since h1 is
  asymmetric, for convolution h1 should be
  reversed, 1 -1).
• So we can then search f(x) for extrema where
• 0 occurs
• -a,b occurs
• a,-b occurs

17
Recall the following from calculus regarding the
second derivative

• Let f be continuous on a,b, and at least twice
  differentiable on (a,b).
• If f(x)gt0 on (a,b), then f is concave up on
  a,b if f(x)lt0 on (a,b), then f is concave
  down on a,b.

18
Recall the second derivative test from calculus

• Let the function f be defined on an open interval
  containing the critical point c where f(c)0,
  and let f be continuous on this interval.
• If f(c)lt0, then c is a local maximum point.
• If f(c)gt0, then c is a local minimum point.
• If f(c)0, then no conclusion is possible
  without further investigation.

19
Also recall

• A point c is called an inflection point of f if
• f is continuous at c
• the graph of f has a tangent line (possibly
  vertical) at c,f(c)
• f is concave up on one side of c and concave down
  on the other side.

20
We can consider the inflection point to be the
location of an edge!
21
Can we calculate the second derivative of an
image?

• Let h1-1, 1.
• Step 1 f(x) f(x) ? h1(x)
• Step 2 f(x) f(x) ? h1(x)
• where ? is cross correlation (convolution).
• So we can then search f(x) for edges where
• 0 occurs
• -a,b occurs
• a,-b occurs

22
Convolution is associative(A ? B) ? C A ?(B ?
C)

• Let h1-1, 1.
• Step 1 f(x) f(x) ? h1(x)
• Step 2 f(x) f(x) ? h1(x)
• where ? is cross correlation (or h11, -1 for
  convolution).
• Given that convolution is associative, is there a
  more computationally efficient way?

23
Can we calculate the second derivative of an
image?

• A better way let h21, -2, 1.
• Then f(x) f(x) ? h2(x)
• Much more computationally efficient.
• But where does h2 come from?
• From -1, 1 ? 1, -1 where ? is cross
  correlation (note the asymmetry), or h1 ? h1
  where ? is convolution.
• See http//en.wikipedia.org/wiki/Cross-correlation
  for differences between convolution and cross
  correlation.

24
Applying masks to images

• Convolution
• Discrete form

Note the reversal of h if h is asymmetric.
25
Frequency response (dB) and phase of first
derivative filter.
26
Frequency response (dB) of first and second
derivative filters.
27
Example from http//homepages.inf.ed.ac.uk/rbf/HIP
R2/zeros.htm
image after zero crossing detection
input image
image after filter
28
Separablefilters

• Danger, Will Robinson!
• Math ahead!

29
Convolution is associative

• (f v) h f (v h), where
• is convolution
• f if the image we wish to convolve
• v is a filter kernel
• h is a another filter kernel

30
Separable

• (f v) h f (v h) (associative)
• If we can express a filter s as
• s v h,
• we can express
• f s f (v h)
• as
• (f v) h

31
Example of separable filter

• The classic Sobel filter
• Convolution of vh is the outer product (or
  tensor product) of two vectors.

32
Outer product of two vectors

• Outer (or tensor) product
• Inner (or dot or scalar) product

33
Outer product of two vectors

• Outer (or tensor) product
• Inner (or dot or scalar) product
• Note The inner product is the trace of the
  outer product.

34
Why do this?

• Consider an MxN image and a PxQ filter kernel.
  The cost of computation is MNPQ.
• For the separable version, the cost is MN(PQ).
• So the savings is PQ/(PQ).
• For a 9x9 kernel, the speed up is about 4.5!

35
Is M separable?

• Note Not all filters are separable!
• How can I tell if a matrix M is separable?
• rank(M) provides an estimate of the number of
  linearly independent rows or columns.
• If the rank(M)1, then it is separable.

36
If separable, how?

• Use svd (singular value decomposition).

37
Octave / Matlab

• s -1 0 1 -2 0 2 -1 0 1 or try ones(5,5) /
  25
• v 1 2 1
• h -1 0 1
• v h result is same as s
• conv2( v, h ) result is same as s
• rank( s ) should be 1
• U,S,V svd( s )
• v2 U(,1) sqrt( S(1,1) ) first col of U
• h2 V(,1)' sqrt( S(1,1) ) first col of VT
• s v2 h2 check - should be 0s
• s v h ditto

38
Octave / Matlab

• How about s 1 2 1 2 4 2 1 2 1 ?
• rank( s ) is 1
• U,S,V svd( s )
• v2 U(,1) sqrt( S(1,1) )
• -1
• -2
• -1
• h2 V(,1)' sqrt( S(1,1) )
• -1 -2 -1
• v2 h2 is s

39
svd singular value decomposition

• M U S VT
• M is an m x n matrix.
• U is an m x m unitary (UTU UUT I) matrix.
• The m cols of U are the left-singular vectors of
  M.
• S is and m x n matrix
• The values along the diagonal of S are the
  singular values of M.
• VT is an n x n matrix
• The n cols of V are the right-singular vectors of
  M.

40
svd and eigenanalysis

• svd is related to eigenanalysis
• The left-singular vectors of M (cols of U) are
  eigenvectors of MMT.
• The right-singular vectors of M (cols of V) are
  eigenvectors of MTM.
• The non-zero singular values of M (diagonal of S)
  are the square roots of the non-zero eigenvalues
  of both MTM and MMT.

41
References

• http//blogs.mathworks.com/steve/2006/10/04/separa
  ble-convolution/
• http//www.gnu.org/software/octave/