% sidex.m - Demonstration of using FFT cross-correlation to compute
% the impulse response of a filter given its input and output.
% This is called "FIR system identification".

Nx = 256;                                   % input signal length
Nh = 100;                                  % filter length
Ny = Nx+Nh-1;                              % max output signal length

% FFT size to accommodate cross-correlation:
Nfft = 2^nextpow2(Nx+Ny-1);                % want power of 2 for FFT

%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 signal
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 of x and y
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 frequency response
plot(1:lge(hxy),real(hxy));

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));