Rizwan Sadiq
International Islamic University, Islamabad Digital Signal Processing LAB
EXPERIMENT # 02: Difference Difference Equations, Convolution, Convolution, Deconvolution Deconvolution and System Identification Identification
Name of Student: Student: …………………………………. …………………………………... Roll No.: …………………………………………… Date of Experiment: ……………………………….. Report submitted on: ……………………………….. ………………………………..
Marks obtained: …………………………………… Remarks: …………………………………………… Instructor’s Signature: ……………………………... Digital Signal Processing Lab
Page 1
Rizwan Sadiq Lab 2 Difference Equation: The difference equation is a formula for computing an output sample at time based on past and present input samples and past output samples in the time domain. We may write the general, causal, LTI difference equation as follows:
FIR System Identification Estimating an impulse response from input-output measurements is called system identification, Cross-correlation can be used to compute the impulse response h(n) of a filter from the crosscorrelation of its input and output signals x(n) and y(n)=x(n)*h(n) , respectively. To see this, note that, by the correlation theorem,
Therefore, the frequency response equals the input-output cross-spectrum divided by the input power spectrum:
where multiplication and division of spectra are defined point wise, i.e.,
.
Digital Signal Processing Lab
Page 2
Rizwan Sadiq %Matlab Code % correlation to compute the impulse response % of a filter given its input and output. % This is called "FIR system identification". Nx = 32; % input signal length Nh = 10; % filter length Ny = Nx+Nh-1; % max output signal length % FFT size to accommodate cross-correlation: Nfft = 2^nextpow2(Ny); % FFT wants power of 2 x = rand(1,Nx); % input signal = noise %x = 1:Nx; % input signal = ramp h = [1:Nh]; % the filter xzp = [x,zeros(1,Nfft-Nx)]; % zero-padded input yzp = filter(h,1,xzp); % apply the filter X = fft(xzp); % input spectrum Y = fft(yzp); % output spectrum Rxx = conj(X) .* X; % energy spectrum of x Rxy = conj(X) .* Y; % cross-energy spectrum Hxy = Rxy ./ Rxx; % should be the freq. response hxy = ifft(Hxy); % should be the imp. response hxy(1:Nh) % print estimated impulse response freqz(hxy,1,Nfft); % plot estimated freq response err = norm(hxy - [h,zeros(1,Nfft-Nh)])/norm(h); disp(sprintf(['Impulse Response Error = ',... '%0.14f%%'],100*err)); err = norm(Hxy-fft([h,zeros(1,Nfft-Nh)]))/norm(h); disp(sprintf(['Frequency Response Error = ',... '%0.14f%%'],100*err));
Convolution in Matlab (without using built in command) Output of LTI system y(n) can be computed by convolving the input signal x(n) with system response h(n). One method to implement is a fast convolution method which includes fft of signal and system response. For much longer convolutions, the savings become enormous compared with ``direct'' convolution. This happens because direct convolution requires on the order of operations (multiplications and additions), while FFT-based convolution requires on the order of N lg( N ) operations, where lg( N ) denotes the logarithm-base-2 of N = 1024;
% FFT much faster at this length
Digital Signal Processing Lab
Page 3
Rizwan Sadiq t = 0:N-1; % [0,1,2,...,N-1] h = exp(-t); % filter impulse reponse H = fft(h); % filter frequency response x = ones(1,N); % input = dc (any signal will do) Nrep = 100; % number of trials to average t0 = clock; % latch the current time for i=1:Nrep, y = conv(x,h); end % Direct convolution t1 = etime(clock,t0)*1000; % elapsed time in msec t0 = clock; for i=1:Nrep, y = ifft(fft(x) .* H); end % FFT convolution t2 = etime(clock,t0)*1000; disp(sprintf([... 'Average direct-convolution time = %0.2f msec\n',... 'Average FFT-convolution time = %0.2f msec\n',... 'Ratio = %0.2f (Direct/FFT)'],... t1/Nrep,t2/Nrep,t1/t2));
Another simple way to implement Convolution A simple function is created using C++ logics to implement convolution formula of sum of products FUNCTION conv(x,h) close all clear all x=input('Enter x: ') h=input('Enter h: ') m=length(x); n=length(h); X=[x,zeros(1,n)]; H=[h,zeros(1,m)]; for i=1:n+m-1 Y(i)=0; for j=1:m if(i-j+1>0) Y(i)=Y(i)+X(j)*H(i-j+1); else end end end Y stem(Y); ylabel('Y[n]'); xlabel('----->n'); title('Convolution of Two Signals without conv function');
Digital Signal Processing Lab
Page 4
Rizwan Sadiq Deconvolution: If the impulse response and the output of the system are known, then the procedure to obtain the unknown input is known as deconvolution. Deconvolution can be implemented by a similar procedure to that of convolution or system identification. i.e y (0) x(0) h(0) And n
y (n) h(m) x (n m ) x(n)
m 1
h(0)
Digital Signal Processing Lab
Page 5
Rizwan Sadiq Lab Tasks
Write code for deconvolution process to find an unknown input signal Write a code for convolution process using convolution as a inner product (hint: use toeplitz matrix or convmtx for help ) In system identification code why we are using zero padding, what will be the output if we don’t use zero padding? Implement a difference equation of IIR filter of any order with coefficients of your own choice with any input signal and impulse response
\ Note: a) This assignment must be submitted before the next lab. b) The assignment submitted must be in proper format as instructed by the teacher to get maximum marks. c) Marks will be deducted on late submissions.
Digital Signal Processing Lab
Page 6
