Download Decimation and Interpolation and more Exercises Advanced Computer Programming in PDF only on Docsity!
Lab#10 Decimation and Interpolation
Downsampling (Decimation):
Decimation is a technique for reducing the number of samples in a discrete-time signal. The operation of
downsampling by factor M describes the process of keeping every Mth sample and discarding the rest.
This is denoted by " ↓ M " in block diagrams.
The decimation process filters the input data with a low pass filter and then resample’s the resulting
smoothed signal at a lower rate.
Decimation is the process of filtering and downsampling a signal to decrease its effective sampling rate,
as illustrated in figure. The filtering is employed to prevent aliasing that might otherwise result from
downsampling.
Why Decimate?
The most immediate reason to decimate is simply to reduce the sampling rate at the output of one system
so a system operating at a lower sampling rate can input the signal. But a much more common motivation
for decimation is to reduce the cost of processing: the calculation and/or memory required to implement a
DSP system generally is proportional to the sampling rate, so the use of a lower sampling rate usually
results in a cheaper implementation.
Code for Decimation:
%%%%%% sequence generation %%%%%%% M=10; n=0:20; x=(1.5).^n; x1=[zeros(1,2000) x];
x2=fliplr(x1); sequence=[x1 x2]; figure(1); subplot(221); stem(sequence); title('sequence graph before decimation') fft_of_sequence=abs(fft(sequence,8192)); x1=[0:1:8191]./(4096); figure(1); subplot(222); plot(x1,fft_of_sequence); title('frequency domain before decimation')
%%%%%%% Prefiltering for decimation %%%%%%%%%% cutoff=1/M; order=400; b=fir1(order,cutoff); % Calculating numirator of FIR filter with cutoff 1/M to avoide alising filtered_sequence=filter(b,1,sequence); % Appling FIR filter to the sequence figure(1); subplot(223); stem(filtered_sequence); title('Sequence graph after filter before decimation') fft_of_filtered_sequence=abs(fft(filtered_sequence,8192)); figure(1); subplot(224); plot(x1,fft_of_filtered_sequence); title('frequency domain after filter before decimation')
%%%%%%% Decimation by factor M %%%%%%%%% deci_sequence=filtered_sequence(1:M:length(filtered_sequence)); % pick sample after gap of M in filterd sequence figure(3); subplot(211); stem(deci_sequence); title('Sequence after filter and decimation') fft_of_deci_sequence=abs(fft(deci_sequence,8192)); figure(3); subplot(212); title('frequency after filter and decimation') plot(x1,fft_of_deci_sequence);
Upsampling (Interpolation):
Interpolation is a technique for increasing the number of samples in a discrete-time signal. The operation
of upsampling by factor L describes the insertion of L-1 zeros between every sample of the input signal.
This is denoted by “↑L “in block diagrams, as in figure.
The process of increasing the sampling rate is called interpolation. Interpolation is upsampling followed
by appropriate filtering. y(n) obtained by interpolating x(n) is generally represented as:
subplot(311); plot(x1,fft_of_sequence); title('Fequency domain before interpolation')
%%%%%%%%%% Up sampling by a factor of L %%%%%%%%% L=2; interpol_sequence =zeros(1,L*(length(sequence))); % Creating zeros vector j = 1;
for i = 1:L:length(interpol_sequence) % adding zeros to the sequence interpol_sequence(1,i) = sequence(1,j); j=j+1; end figure(1); subplot(312),stem(interpol_sequence);title('Sequence after interpolation') fft_of_interpol_sequence = abs(fft(interpol_sequence,131072)); figure(2); subplot(312); plot(x1,fft_of_interpol_sequence);title('fequency domain after interpolation')
MATLAB COMMANDS:
1) Decimate:
-It reduces the sample rate.
y = decimate(x, r)
y = decimate(x,r,n)
y = decimate(x,r,'fir')
- y = decimate(x, r) reduces the sample rate of x by a factor r.
- The decimated vector y is r times shorter in length than the input vector x.
- By default, decimate employs an eighth-order lowpass Chebyshev Type I filter.
- It filters the input sequence in both the forward and reverse directions to remove all phase distortion,
effectively doubling the filter order.
- y = decimate(x,r,n) uses an order n Chebyshev filter. Orders above 13 are not recommended because of
numerical instability.
- y = decimate(x,r,'fir') uses a 30-point FIR filter, instead of the Chebyshev IIR filter.
- Here decimate filters the input sequence in only one direction. This technique conserves memory and is
useful for working with long sequences.
2) Interp:
-It increases the sampling rate.
y = interp(x, r)
y = interp(x, r, l, alpha)
[y, b] = interp(x, r, l, alpha)
- y = interp(x, r) increases the sampling rate of x by a factor of r. The interpolated vector y is r times
longer than the original input x.
- y = interp(x, r, l, alpha) specifies l (filter length) and alpha (cut-off frequency). The default value for l is
4 and the default value for alpha is 0.5.
- [y,b] = interp(x,r,l,alpha) returns vector b containing the filter coefficients used for the interpolation.
3) Fliplr:
B = fliplr (A)
-It returns A with columns flipped in the left-right direction, that is, about a vertical axis.
- If A is a row vector, then fliplr (A) returns a vector of the same length with the order of its elements
reversed.
- If A is a column vector, then fliplr (A) simply returns A.
Lab Tasks
Task-1:
a) Generate a sinusoid of frequency f0 = 3kHz and sampling frequency fs=8kHz using
x[n] = sin((2 * pi * f0 / fs))*n) where n = 0:1: 15999.
Stem the Fourier transform of x[n] using fft command (zero-frequency component should be at center.
b) Now decimate it by M = 2 & if you see the Fourier transform of this sequence, it will show a sinusoid
at 1kHz which is an aliased component.
c) Now pass x[n] through an anti-aliasing filter of wc = pi/2 and see the spectrum. You will see that the
signal contains no frequency components as wc = pi/2 =
2kHz and its decimated signal will also contain no frequency components. Hence aliasing is avoided with
the help of anti-aliasing filter.
Task-2:
Apply all three steps given in task-1 to following sequences:
1) x1[n] = cos((2 * pi * f0 / fs))n) + cos((2 * pi * f1 / fs))n)
where n = 0:1: 15999, f0 = 2kHz, f1 = 6kHz, fs = 10kHz
2) x2[n] = cos((2 * pi * f0 / fs))n) + cos((2 * pi * f1 / fs))n) + cos((2 * pi * f2 / fs))*n)
where n = 0:1: 15999, f0 = 2kHz, f1 = 4kHz, f1 = 3kHz, fs = 10kHz
