Transcript Chapter 12

Chapter 12
Case Studies in Digital Signal
Processing
Case Studies
• Dual-Tone Multifrequency (DTMF)
Signaling
• A Software Impedance Bridge
• Details of JPEG Compression
• FBI Fingerprint Image Compression Using
Wavelets
Case I: Dual-Tone Multifrequency Signaling
• Used to communicate Touchtone® pad
information via telephone links
• Based on CCITT Q.23 and Q.24 standards
DTMF Frequency Set
4 Low Frequencies
3 High Frequencies
No harmonic
relationships
Example: “5” is the
sum of 770 Hz and
1336 Hz
Operational Specifications for
DTMF Signaling
Parameter
Operation
Specification
Frequency Tolerance
Accept
< = 1.5%
Reject
>= 3.5%
Accept
> = 40 ms
Reject
< 23 ms
Signal Interruption
Reject as two separate tones
< 10 ms
Twist
Forward - Accept
<= 8 dB
Reverse - Accept
<= 4 dB
Signal Duration
Tone Generation
bL  sin( L )
bH  sin( H )
aL  2 cos( L )
aH  2 cos( H )
H
(bL  bH ) z 1  (aH bL  aLbH ) z 2  (bL  bH ) z 3
1  (aL  aH ) z 1  (aL aH  2) z 2  (aL  aH ) z 3  z 4
(bL  bH ) z 1  (aH bL  aLbH ) z 2  (bL  bH ) z 3
H
1  (aL  aH ) z 1  (aL aH  2) z 2  (aL  aH ) z 3  z 4
The transfer function is a 4-pole resonator with
frequencies at one set of the tone pairs (ΩL and ΩH)
Example: Digit “7”
(852 and 1209 Hz)
>> f1=852;
>> f2=1209;
>> fs=8000;
>> bL=sin(2*pi*f1/fs);
>> bH=sin(2*pi*f2/fs);
>> aL=2*cos(2*pi*f1/fs);
>> aH=2*cos(2*pi*f2/fs);
>> b=[0,bL+bH,-(aL*bH+aH*bL),bL+bH];
>> a=[1,-(aL+aH),aL*aH+2,-(aL+aH),1];
>> x=[1,zeros(1,2000)];
%This is an impulse signal to
start the oscillator
>> y=filter(b,a,x);
>> dtft_demof(y,0,1400,1024,8000); %The sampling
frequency is 8 kHz
>> figure, zplane(b,a) % Examine the poles and zeros
Discrete Time Fourier Transform
1000
1
900
0.8
800
0.6
700
Imaginary Part
0.4
600
500
400
0.2
0
-0.2
-0.4
300
-0.6
200
-0.8
100
-1
0
0
200
400
600
800
Hz
1000
1200
1400
-1
-0.5
0
Real Part
0.5
1
Tone Detection and Validation
•
•
•
•
•
Rejection of dial-tone
Estimation of spectral content (Goertzel
Algorithm)
Estimation of second-harmonic power
Determination of forward and reverse twist
ratios
Tone duration check
Rejection of Dial Tone
(350 Hz and 440 Hz )
Cascaded Notch Filter
>> fs=8000;
% The system
sampling frequency is 8 Hz
>> f350=2*pi*350/fs;
% Compute the
digital frequencies for each notch
>> f440=2*pi*440/fs;
>> deltaf=2*pi*5/fs;
% Let the -3 dB
width of each notch be 5 Hz
>> r=1-deltaf/2;
% Compute the pole
radius
>> g350=abs(1-2*r*cos(f350)+r^2)/(2*abs(1cos(f350))); % Compute the “g” factors
>> g440=abs(1-2*r*cos(f440)+r^2)/(2*abs(1cos(f440)));
>> b350=g350*[1,-2*cos(f350),1]; % Compute
the b and a coefficients for each filter
>> b440=g440*[1,-2*cos(f440),1];
>> a350=[1,-2*r*cos(f350),r^2];
>> a440=[1,-2*r*cos(f440),r^2];
>> b2=conv(b350,b440);
% Convolve
the coefficient vectors
>> a2=conv(a350,a440);
>> x=analog([350,440,697,1209],[1,1,1,1],1000,8000);
>> y=filter(b2,a2,x);
>> subplot(2,1,1),dtft_demof(x,0,1400,2048,8000);
title('Digit 1 with Dial-Tone')
>> subplot(2,1,2),dtft_demof(y,0,1400,2048,8000);
title('Digit 1 with Dial-Tone Removed')
Digit 1 with Dial-Tone
4000
3000
2000
1000
0
0
200
400
600
800
1000
Hz
Digit 1 with Dial-Tone Removed
0
200
400
600
1200
1400
1200
1400
4000
3000
2000
1000
0
800
Hz
1000
Determination of Tone Frequencies
• Tone frequencies are determined by the Goertzel
algorithm
• The Goertzel algorithm is based on computing the DFT
at just the expected tone frequencies and determining
which frequencies show the largest response
• It can be shown that the Goertzel algorithm is more
efficient than the FFT if M < log2N where M is the
number of frequencies to be computed by the DFT and
N is the signal length.
• If M = 7, then N = 136, then the Goertzel algorithm can
operate on 136 samples of the MF tone signal,
• At a sampling frequency of 8 kHz, tone detection
requires a minimum of 136*.125 = 17 ms.
Goertzel Algorithm Example
(20 Sample Signal, 250 Hz and 500 Hz)
N- 1
X[ k] =
å
x[n]e-
j 2π( k / N ) n
,
k = 0,1, 2,
n= 0
,N- 1
W1 = 2π
k1
3
= 2π
N
20
W2 = 2π
k2
5
= 2π
N
20
X[k] = X(Wk )
Wk = 2π
W= 2π
k
,
N
k = 0 ,1, 2 ,
f
k
= 2π
fs
N
fN
fs
( 250)( 20)
k1 =
= 2.5 @ 3
2000
(500)( 20)
k2 =
=5
2000
k=
,N- 1
19
X [k1 ]   x[n]e
2
3
n
20
2
5
n
20
n 0
19
X [k2 ]   x[n]e
n 0
Goertzel Algorithm Example
(20 Sample Signal, 250 Hz and 500 Hz)
>> n=1:20;
>> f1=2*pi*250/2000;
>> f2=2*pi*500/2000;
>> x=sin(f1*n)+sin(f2*n);
>> dtft_demo(x,0,2*pi,512); % Custom M-file for the
DTFT
>> hold
>> [X,omegax]=dft_demo(x); % Custom M-file for the
DFT
>> stem(omegax/pi,abs(X), 'k');
>> title('250 and 500 Hz, DTFT (Curve) and DFT (Stem)')
>> hold off
>> dtft_demo(x,0,2*pi,512);
>> hold
>>
stem([omegax(4),omegax(6)]/pi,[abs(X(4)),abs(X(6)
)],'k')
>> title('DTFT (Curve) and Goertzel Results (Stem)')
>> hold off
DTFT (Curve) and Goertzel Results (Stem)
250 and 500 Hz, DTFT (Curve) and DFT (Stem)
12
12
10
10
8
8
6
6
4
4
2
2
0
0
0.2
0.4
0.6
0.8
1
1.2
Units of Pi
1.4
1.6
1.8
2
0
0
0.2
0.4
0.6
0.8
1
1.2
Units of Pi
1.4
1.6
1.8
2
“goertzel_vector” M-file
>> help goertzel_vector
G=GOERTZEL_VECTOR(N,FS,F)
This function produces a vector of index values, G, to use in the Goertzel
algorithm. The index values correspond to the closest Hertzian frequencies
given
in the vector F, while N = number of sample values in the input signal to
the Goertzel command and FS is the sampling frequency of the input signal.
>> f=[697,770,852,941,1209,1336,1477];
>> g=goertzel_vector(136,8000,f)
g=
13 14 15 17 22 24 26
Goertzel index set for the
DTMF tones when N = 136
and Fs = 8 kHz
Goertzel Algorithm Example
Detection of Digit “1”
>> f=[697,770,852,941,1209,1336,1477];
>> g=goertzel_vector(136,8000,f);
>> n=0:135;
>> f1=2*pi*697/8000;
>> f2=2*pi*1209/8000;
>> tone1=sin(f1*n)+sin(f2*n);
>> X1=goertzel(tone1,g);
>> abs(X1);
>> stem(f,abs(X1))
>> xlabel('Hz')
>> title('Goertzel Algorithm Output for the Tone
"1" (697 and 1209 Hz)')
Goertzel Algorithm Output for the Tone "1" (697 and 1209 Hz)
70
60
Each stem corresponds to
one of the 7 tone
frequencies. The
maximum response
decodes the dual tone as
697 Hz and 1209 Hz = digit
“1”
50
40
30
20
10
0
600
700
800
900
1000
1100
Hz
1200
1300
1400
1500
Goertzel Algorithm Example
Rejection of Digit “0”
>> n=0:135;
>> f1off=2*pi*1001/8000;
>> f2off=2*pi*1396/8000;
>> tone0off=sin(f1off*n)+sin(f2off*n);
>> X0off=goertzel(tone0off,g);
>> stem(f,abs(X0off))
>> xlabel('Hz')
>> title('Goertzel Output for "0" Tones Offset by
+60 Hz')
Goertzel Output for "0" Tones Offset by +60 Hz
25
The valid tones for “0” are
941 and 1336 Hz. Here the
tones have been offset by
+60 Hz (3.5%). The upper
tone is detected but no
lower tone is detected
greater than 4 dB below the
upper tone. Result:
rejection as a valid tone set.
20
15
10
5
0
600
700
800
900
1000
1100
Hz
1200
1300
1400
1500
Case II: Identification of Physical Impedances
Using an Adaptive Filter
+
Unknown System H(z)
ADC
x(t) = sin(? t)
ADC
x[k]
y[k]
S
-
ek
Wiener Filter
Wiener Filter

N 1
y[k ]   w[i ]x[k  i ]
i 0
Estimate of
unknown system
The unknown system changes the amplitude and phase of
an input sinusoid. The Wiener filter models the response
as an FIR system with transfer function H(z). From the
frequency response of H(z), the unknown impedance can
be estimated.
Analog Impedance Circuit
®
Vx
®
Vr
®
Zx
=
®
Rm + Z x
=
Rx + jX
Rm + Rx + jX
®
æ Rx + jX ö÷
ç
÷= F ( Z ) Vr = AÐθ
Vx = Vr ç
÷
çè Rm + Rx + jX ÷
ø
®
®
The unknown impedance will change
the amplitude and phase of an input
sinusoid by A and θ respectively. The
resistance Rm is assumed to be known.
F(Z) is the “impedance factor” that
results in A and θ.
Wiener Filter Estimate of H(z) and the Unknown
Impedance
In the z-domain:
Vx ( z) = H( z)Vr ( z)
H(W) = H( z) z= e jW
W= 2π
Since both H(Ω) and the impedance factor
F(Z) are complex numbers with the same
magnitude and phase angle, we can
conclude:
®
F (Z ) = H (W) =
f
fs
Rx + jX
Zx
=
Rm + Rx + jX R + Z®
m
x
Using the Wiener filter coefficients:
A = H(W)
θ = ÐH(W)
H (W) = w0 + w1e®
jW
+ w2e-
j 2W
+
Rm H (W)
Zx =
= Rx + jX
1- H (W)
+ wN- 1e-
j ( N - 1)W
Example
>> Rm=100;
>> Rx=50;
% This is the “unknown” impedance
>> X=-125;
>> FZ=(Rx+j*X)/(Rm+Rx+j*X); % This is the complex response of the physical
% impedances
>> amp=abs(FZ)
% This gives the amplitude response
amp =
0.6895
>> deg=angle(FZ)*180/pi % This gives the phase response in degrees
deg =
-28.3930
>> n=0:44099;
>> omega=2*pi*200/44100;
%The digital frequency for 200 Hz at the
%sampling rate of 44.1 kHz.
>> x=sin(omega*n);
% This is 1 second of the sampled excitation signal
>> y=amp*sin(omega*n+angle(FZ));
%This is the digital output signal with an
% amplitude and phase that would be
% produced by the unknown system
% response
Example, Continued
>> [err,n_hat,w]=adapt2(6,x,y);
%Run the adaptive filter routine
>> plot(err); title('Error vs. Samples'); xlabel('Sample')
>> b=w';
%Put the estimated impulse response into a row vector
>> z=exp(j*omega); %The complex unit circle value of the input frequency
>> H=H_evaluate(b,1,z); %Use the custom M-file H_evaluate to compute the
% response of the system at the input frequency
>> Rm*H/(1-H)
% Compute the real and imaginary parts of the unknown
% impedance
ans =
5.0002e+001 -1.2500e+002
% The estimate of the unknown impedance Rx +jX
Example, Continued
Error vs. Samples
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4
0
0.5
1
1.5
2
2.5
Sample
3
3.5
4
4.5
4
x 10
The error signal shows the convergence
of the Wiener filter to the estimate of the
unknown system.
Case III: JPEG Compression
• Joint Photographic Experts Group standard
adopted in 1994 as ISO 10918-1 (International
Organization for Standards).
• JPEG is an example of transform compression
using the Discrete Cosine Transform.
• Because some of the transform information is
deleted in the algorithm, JPEG is termed lossy
compression.
• Suitable for photographs but less so for line
drawings and textual graphics.
The Concept of Transform Compression
The Fourier series
approximations to a
triangle wave only
require a few low
frequency components
for an adequate
representation. The 10term signal is a 10-to-1
“compression” of the
100-term signal.
Discrete Cosine Transform
• The DCT is a variation on the discrete Fourier
Transform (DFT)
• Provides a real-valued frequency
representation in the range Ω = [0,π].
  (2n  1)(k  1) 
y[k ]  w[k ] x[n]cos 
 .......k  1,..., N
2
N


n 1
1
w[k ] 
...k  1
N
N
w[k ] 
2
...2  k  N
N
Comparison of DFT and DCT for a
Ω = π/4 Sinusoid
>> n=0:20;
>> f=pi/4;
>> x=sin(f*n);
>> yf=fft(x);
>> yc=dct(x);
>> omegaf=0:2*pi/20:2*pi;
>> omegac=0:pi/20:pi;
>> subplot(1,2,1),stem(omegaf/pi,abs(yf)),title('Discrete
Fourier Transform')
>> xlabel('units of pi')
>> subplot(1,2,2),stem(omegac/pi,yc),title('Discrete
Cosine Transform')
>> xlabel('units of pi')
Discrete Fourier Transform
Discrete Cosine Transform
8
2
7
1.5
1
6
The DCT has only real
coefficients and is defined only
in the interval [0,π] and is
therefore a more “efficient”
transform of the sinusoid.
0.5
5
0
4
-0.5
3
-1
2
-1.5
1
0
-2
0
0.5
1
1.5
units of pi
2
-2.5
0
0.5
units of pi
1
Two-Dimensional DCT
  (2m  1) p 
  (2n  1) q 
Y [ p, q]  u[ p]v[q]   x[m, n]cos 
cos



2M
2N




m 0 n 0
p  0,1, 2,...M  1
M 1 N 1
q  0,1, 2,...N  1
1
u[ p ] 
... p  0
M
2
....1  p  M  1
M
1
v[q] 
...q  0
N
u[ p ] 
v[q] 
2
...1  q  N  1
N
Steps in JPEG Compression
•
•
•
•
•
•
Convert the image matrix pixel values to the correct numerical storage
class for arithmetic computations. After this step the pixel values (matrix
values) will range from 0 to 1 in grayscale.
Level-shift the matrix values so that the middle grayscale value is zero.
This would mean that 0.5 is subtracted from each pixel value. This
basically eliminates any “DC bias” in the values and makes the DCT
more sensitive to frequency variations.
Compute the 2-dimensional DCT on the each 8-by-8 block of the levelshifted image matrix.
Multiply the DCT matrix block by an 8-by-8 “mask” matrix, element-byelement, to select a variable number of DCT coefficients, starting with the
lowest frequency coefficients. This is the step that actually produces the
compression of the image.
Convert the masked DCT block to a vector by means of a “zigzag”
reading of the matrix values and discard trailing zeros. The block vectors
form the compressed file for the image.
Reconstruct the image by reversing the above steps. The key step is
taking the inverse DCT of the masked coefficient blocks.
Step 1 – Create the 8-by-8 block and convert to
double precision
>> I=imread('pout.tif');
>> size(I)
ans =
291 240
>> block=I(101:108,101:108)
block =
120 123 134 142 144 146 149 146
119 123 136 144 147 144 148 145
117 122 133 139 147 144 147 143
117 121 128 133 144 144 146 144
112 119 127 130 138 142 145 145
106 112 121 123 128 133 139 142
95 101 109 112 114 117 125 130
91 91 102 107 108 108 108 113
>> b2=im2double(block) % Convert the matrix values to storage class “double”
b2 =
0.4706 0.4824 0.5255 0.5569 0.5647 0.5725 0.5843 0.5725
0.4667 0.4824 0.5333 0.5647 0.5765 0.5647 0.5804 0.5686
0.4588 0.4784 0.5216 0.5451 0.5765 0.5647 0.5765 0.5608
0.4588 0.4745 0.5020 0.5216 0.5647 0.5647 0.5725 0.5647
0.4392 0.4667 0.4980 0.5098 0.5412 0.5569 0.5686 0.5686
0.4157 0.4392 0.4745 0.4824 0.5020 0.5216 0.5451 0.5569
0.3725 0.3961 0.4275 0.4392 0.4471 0.4588 0.4902 0.5098
0.3569 0.3569 0.4000 0.4196 0.4235 0.4235 0.4235 0.4431
Step 2 – Level-shift the matrix values
>> b3=b2-.5
b3 =
-0.0294
-0.0333
-0.0412
-0.0412
-0.0608
-0.0843
-0.1275
-0.1431
% Subtract 0.5 from each matrix value to create +/- 1 pixel values
-0.0176
-0.0176
-0.0216
-0.0255
-0.0333
-0.0608
-0.1039
-0.1431
0.0255 0.0569 0.0647 0.0725
0.0333 0.0647 0.0765 0.0647
0.0216 0.0451 0.0765 0.0647
0.0020 0.0216 0.0647 0.0647
-0.0020 0.0098 0.0412 0.0569
-0.0255 -0.0176 0.0020 0.0216
-0.0725 -0.0608 -0.0529 -0.0412
-0.1000 -0.0804 -0.0765 -0.0765
0.0843
0.0804
0.0765
0.0725
0.0686
0.0451
-0.0098
-0.0765
0.0725
0.0686
0.0608
0.0647
0.0686
0.0569
0.0098
-0.0569
Step 3 – Compute the 2-D DCT of the level-shifted
matrix
>> c1=dct2(b3)
c1 =
0.0059
0.3411
-0.1531
0.0645
-0.0186
-0.0089
0.0019
-0.0031
-0.3088
-0.0007
0.0294
-0.0277
0.0098
-0.0089
0.0027
-0.0031
-0.0867
-0.0428
-0.0058
0.0242
-0.0030
0.0143
-0.0042
-0.0009
-0.0315
0.0175
-0.0161
-0.0099
0.0160
0.0024
0.0051
-0.0037
-0.0010
-0.0007
0.0146
-0.0132
0.0020
-0.0024
-0.0007
-0.0046
0.0087
0.0093
0.0154
0.0017
-0.0026
-0.0025
0.0004
-0.0005
0.0009
-0.0100
0.0056
0.0003
0.0055
-0.0042
0.0039
-0.0019
0.0203
0.0059
-0.0076
-0.0064
-0.0037
0.0001
0.0005
-0.0006
Step 4 – Select only low index values of the
DCT matrix with a mask.
>> mask=mask8(4)
mask =
1 1
1 1
1 1
1 0
0 0
0 0
0 0
0 0
1
1
0
0
0
0
0
0
>> c2=c1.*mask
1
0
0
0
0
0
0
0
% This is a custom M-file for mask generation.
% It will select a triangle-shaped portion of the low index
% values of the DCT matrix.
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
% Multiply the DCT matrix by the mask element-by-element
c2 =
0.0059 -0.3088 -0.0867 -0.0315
0.3411 -0.0007 -0.0428
0
-0.1531 0.0294
0
0
0.0645
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Step 5 – Read out the retained matrix
values into a vector using a “zig-zag”
>> c2
c2 =
0.0059 -0.3088 -0.0867 -0.0315
0.3411 -0.0007 -0.0428
0
-0.1531 0.0294
0
0
0.0645
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
>> c2vector=zz(c2);
>> length(c2vector)
ans =
64
>> c2vector(1:10)'
ans =
0.0059
-0.3088
0.3411
-0.1531
-0.0007
-0.0867
-0.0315
-0.0428
0.0294
0.0645
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
The vector of the
retained matrix
represents the
“compressed”
information for the
.jpg file
Step 6a – Reconstruct the matrix by computing the
inverse DCT of the masked coefficient matrix
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
•
>> c2=dezz(c2vector) % This command “de-zig-zags” the stored zig-zag vector
ans =
0.0059 -0.3088 -0.0867 -0.0315
0
0
0
0
0.3411 -0.0007 -0.0428
0
0
0
0
0
-0.1531 0.0294
0
0
0
0
0
0
0.0645
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
>> r1=idct2(c2)
r1 =
-0.0312 -0.0043 0.0331 0.0619 0.0745 0.0754 0.0733 0.0722
-0.0396 -0.0130 0.0242 0.0536 0.0677 0.0709 0.0711 0.0715
-0.0473 -0.0214 0.0151 0.0449 0.0612 0.0680 0.0720 0.0748
-0.0498 -0.0255 0.0092 0.0384 0.0561 0.0665 0.0745 0.0800
-0.0567 -0.0347 -0.0032 0.0237 0.0414 0.0540 0.0652 0.0729
-0.0807 -0.0614 -0.0340 -0.0106 0.0055 0.0187 0.0317 0.0409
-0.1187 -0.1018 -0.0783 -0.0586 -0.0446 -0.0319 -0.0182 -0.0082
-0.1487 -0.1332 -0.1120 -0.0946 -0.0822 -0.0701 -0.0562 -0.0460
Step 6b – Continue reconstruction by restoring the
pixel levels
>> r2=r1+.5
r2 =
0.4688
0.4604
0.4527
0.4502
0.4433
0.4193
0.3813
0.3513
0.4957
0.4870
0.4786
0.4745
0.4653
0.4386
0.3982
0.3668
0.5331
0.5242
0.5151
0.5092
0.4968
0.4660
0.4217
0.3880
0.5619
0.5536
0.5449
0.5384
0.5237
0.4894
0.4414
0.4054
0.5745
0.5677
0.5612
0.5561
0.5414
0.5055
0.4554
0.4178
0.5754
0.5709
0.5680
0.5665
0.5540
0.5187
0.4681
0.4299
0.5733
0.5711
0.5720
0.5745
0.5652
0.5317
0.4818
0.4438
0.5722
0.5715
0.5748
0.5800
0.5729
0.5409
0.4918
0.4540
Step 6c – Compare to the original image block
matrix
>> b2-r2
ans =
0.0018
0.0063
0.0061
0.0086
-0.0041
-0.0036
-0.0087
0.0056
-0.0133 -0.0076 -0.0051 -0.0098 -0.0029 0.0110
-0.0046 0.0091 0.0111 0.0088 -0.0062 0.0093
-0.0001 0.0065 0.0002 0.0153 -0.0033 0.0045
-0.0000 -0.0072 -0.0168 0.0086 -0.0018 -0.0019
0.0013 0.0013 -0.0139 -0.0002 0.0029 0.0035
0.0006 0.0085 -0.0070 -0.0036 0.0028 0.0134
-0.0021 0.0057 -0.0022 -0.0083 -0.0093 0.0084
-0.0099 0.0120 0.0142 0.0058 -0.0064 -0.0202
0.0004
-0.0029
-0.0140
-0.0152
-0.0043
0.0159
0.0180
-0.0108
Analysis - Compute the average percent error between the original and reconstructed
block:
>> mean2(((b2-r2)./b2)*100)
ans =
-0.0314
Full Image JPEG:
The Block Processing Command
>> help blkproc
BLKPROC Implement distinct block processing for image.
B = BLKPROC(A,[M N],FUN) processes the image A by applying the function
FUN to each distinct M-by-N block of A, padding A with zeros if
necessary. FUN is a function that accepts an M-by-N matrix, X, and
returns a matrix, vector, or scalar Y:
Y = FUN(X)
BLKPROC does not require that Y be the same size as X. However, B is
the same size as A only if Y is the same size as X.
B = BLKPROC(A,[M N],FUN,P1,P2,...) passes the additional parameters
P1,P2,..., to FUN.
Full Image JPEG:
The Compression M-File
>> type cjpeg1
function cx = cjpeg1(x,mask)
% function cx = cjpeg1(x,mask)
%
% Basic jpeg-like transform coder.
% MASK is an 8 by 8 matrix that sets the DCT coefficients
% to be saved (i.e., compression amount), generated by the M-file MASK8.M.
% <x> is the input image.
% The retained coefficients are returned in <cx>.
% The DC coefficient should be top left of each 8x8 block
% This command uses DCT2 for computing the block transform.
x=im2double(x); %convert image to double precision for the inline commands
x = x-.5;
%level-adjust the image around the mid-range grayscale
fun=@dct2;
%Define the 2-D DCT function
B = blkproc(x,[8 8],fun); %Block process the image to create DCT blocks
cx=blkproc(B, [8 8], 'x.*P1',mask); %Multiply the blocks by the mask
%end of cjpeg1
Full Image JPEG:
M-Files
cjpeg1 – This function carries out an 8-by-8 sub-block computation of
the DCT on a level-adjusted black-and-white image and then
multiplies each block by an 8-by-8 mask generated by mask8.
mask8 – This file generates an 8-by-8 logical mask matrix that will
select the DCT coefficients from the lowest to the highest index
values. When multiplied by an 8-by-8 DCT block, this selects just
a subset of the coefficients in the block, setting the others to zero.
This effects the compression of the block.
djpeg1 – This function reconstructs (decompresses) an image from the
compressed transform file generated by cjpeg.
jpeg_demo – This combines the compression and reconstruction
routines (cjpeg1 and djpeg1) and generates an absolute-value
error image between the original image and its reconstruction. It
also computes the compression ratio by determining the fraction
of DCT coefficients retained after the block masking operation.
Full Image JPEG:
Example
>> I=imread('moon.tif');
>> mask=mask8(4); % 6-to-1
compression
>> jpeg_demo(I,mask);
Original Image
Compression 6 to 1
6 to 1 Compression
>> mask=mask8(2); % 21-to-1
compression
>> figure, jpeg_demo(I,mask);
10x Absolute Error
Original Image
Compression 21 to 1
10x Absolute Error
21 to 1 Compression
Case IV: FBI Fingerprint Compression Using
Wavelets
• FBI fingerprint information could require up to
300 terabytes of storage.
• Compression of the image data by factors up to
20-to-1 were desired.
• JPEG compression is not satisfactory because
of the loss of fine image detail.
• Wavelet compression gives superior results at
high compression ratios.
• The method used is similar to the new JPEG
2000 image compression standard.
Uncompressed Fingerprint Image
>> load detfingr.mat
>> FPbw=X/(max(max(X)));
>> imshow(FPbw),title('Uncompressed B/W Intensity Fingerprint Image')
Uncompressed B/W Intensity Fingerprint Image
JPEG Compressed Fingerprint Images
>> mask_2=mask8(2); %Set a DCT mask to 21-to-1 compression
>> mask_3=mask8(3); %Set a DCT mask to 11-to-1 compression
>> FPr_11=jpeg_demo(FPbw,mask_3); %Generate a reconstructed image 11-tto-1
>> close %Close the figure created by jpeg_demo
>> FPr_21=jpeg_demo(FPbw,mask_2); % Generate a reconstructed image 21-to-1
>> close %Close the figure created by jpeg_demo
>> imshow(FPr_11),title('Reconstructed JPEG 11-to-1 Compression')
>> figure,imshow(FPr_21),title('Reconstructed JPEG 21-to-1 Compression')
Reconstructed JPEG 11-to-1 Compression
Reconstructed JPEG 21-to-1 Compression
JPEG Compressed Fingerprint Images
(21-to-1 Compression)
Reconstructed JPEG 21-to-1 Compression
Note the loss of fine
detail and the block
artifacts in the
reconstructed image
Wavelet Compression Using the “wdencmp “
Command (help file)
[XC,CXC,LXC,PERF0,PERFL2]=WDENCMP('gbl',X,'wname',N,THR,SORH,KEEPAPP)
returns a de-noised or compressed version XC of input
signal X (1-D or 2-D) obtained by wavelet coefficients
thresholding using global positive threshold THR.
Additional output arguments [CXC,LXC] are the
wavelet decomposition structure of XC,
PERFL2 and PERF0 are L^2 recovery and compression
scores in percentages.
PERFL2 = 100*(vector-norm of CXC/vector-norm of C)^2
where [C,L] denotes the wavelet decomposition structure
of X.
Wavelet decomposition is performed at level N and
'wname' is a string containing the wavelet name.
SORH ('s' or 'h') is for soft or hard thresholding
(see WTHRESH for more details).
If KEEPAPP = 1, approximation coefficients cannot be
thresholded, otherwise it is possible.
Wavelet Compression Using the “wdencmp”
>> wav='bior4.4'; %Set the wavelet family to be used
>> level=3;
% Set the decomposition level
>> sorh='h';
% Set the thresholding method to “hard”
>> keepapp=1;
% Do not threshold the approximation coefficients
>> thr=40;
% Set the thresholding level and use globally
(‘gbl’)
>> [XC,CXC,LXC,PERF0,PERFL2]=
wdencmp('gbl',X,wav,level,thr,sorh,keepapp);
>> imshow(XC/max(max(XC))),colormap(bone(128));
>> title('20-to-1 Compression, bior4.4 Wavelet, Threshold = 40')
>> PERF0
PERF0 =
94.8699
Wavelet Compression Using wdencmp
Results for bior4.4 Wavelet Family
20-to-1 Compression, bior4.4 Wavelet, Threshold = 40
Wavelet Compression Using jpeg2000_demo
Results for haar Wavelet Family
>> load detfingr.mat; % Load the image
>> wav='haar';
% Set the wavelet family
>> level=3;
% Set the decomposition level
>> thresh=40;
% Set the thresholding for compression
>> [XC,cmp_ratio,energy,thr]=jpeg2000_demo(X,wav,level,thresh);
>> close
% Close the figure created by jpeg2000_demo
>> imshow(XC/max(max(XC))),colormap(bone(128)); %Display
intensity image
>> title('19-to-1 Compression, Haar Wavelet, Threshold =40')
>> cmp_ratio
% Compression ratio is 19-to-1 for this example
cmp_ratio =
19
Wavelet Compression Using jpeg2000_demo
Results for haar Wavelet Family
19-to-1 Compression, Haar Wavelet, Threshold =40
Summary
• Dual Tone Multifrequency Signaling
– Method by which tone-coded information is exchanged in the telephone
network
– Digital oscillators used for tone generation
– The Goertzel algorithm used for tone identification
• DSP Impedance Bridge
– Adaptive filter system identification is used to model the magnitude and
phase response of the Thevenin equivalent of a physical impedance.
• JPEG Compression
– A form of block-by-block transform compression of images using the
discrete cosine transform (DCT)
• Fingerprint Image Compression
– The use of wavelets for image compression
– Similar to the new JPEG 2000 image compression standard