%
%  Transfer Function Estimate of a Linear System  
%
Fs = 1000;                             % sampling frequency
T = 1;                                 % one second worth of samples�@�@�@�@�@
t = (0:1/Fs:T);                        % time vector
x = sin(2*pi*50*t) + sin(2*pi*120*t) ; % + 0.1*randn(size(t));    
% x = square(2*pi*4*t);                % square wave  
subplot(3,2,1); plot(t,x)              %                  input signal 
axis([0 T/4 -2 2]); legend('input signal x')
%
%                                        < Welch's method of PSD estimation>
nfft = 256;                                     % length of FFT
window = hanning(256);                          % window function
noverlap = 128;                                 % number of samples overlap
dflag = 'none';                                 % detrending option
[Pxx,f] = psd(x,nfft,Fs,window,noverlap,dflag); % Power Spectral Density 
subplot(3,2,3); plot(f,10*log10(Pxx))           % PSD of input signal (dB)
legend('Pxx')
%
h = ones(1,10)/10;                    % moving average filter
Ht = freqz(h,1,f,Fs);                 % theoretical frequency response
subplot(3,2,5); plot(f,abs(Ht)); 
legend('theoretical response')
%
y = filter(h,1,x);                    % filtering: y=conv(h,x);
subplot(3,2,2); plot(t,x)             %                  output signal 
axis([0 T/4 -2 2]); legend('output signal y')
[Pxy,f] = csd(x,y,nfft,Fs,window,noverlap,dflag);  % Cross Spectral Density
[Pyy,f] = psd(y,nfft,Fs,window,noverlap,dflag);
subplot(3,2,4); plot(f,10*log10(abs(Pyy)))         % PSD of output signal (dB)
legend('Pyy')
%
Hc = Pxy./Pxx;                        % estimated frequency response                             
subplot(3,2,6); plot(f,abs(Hc)); 
legend('estimated response')       
% ---------------------------------------------------------------------------------------
%                                                                2001.7.02  by K.Tsukada