% demonstrating the error in approximating the variance of
% the periodogram for a more general wide sense stationary process


L=10; % filter length
N=10000; % data length
sigma=1;
h=rand(1,L);  % generating a FIR filter with random pos. components

e= randn(1,N)*sigma; % gen. input white noise with power sigma^2

x=conv(h,e);  % gen. filter output - conv(e,h) gives same answer
x=x(1:N); % discard last transient part

xwin=x(ceil(N/2)+1:N); % take second half as windowed output

ewin=e(ceil(N/2)+1:N-L+1);  % take same part of input signal; 
                            %take L-1 less data because convolution will increase the lenght by L-1

xwinApprox=conv(h,ewin);  % convolve filter with windowed input
%xwinApprox= xwinApprox(1:length(ewin)); % throw away last transient part

% compute power spectrums using the periodogram
Ndata=length(xwin);
res=4*Ndata;  % dft resolution (implemented with zero padding), ensure res >= N
pwrSpecXwin=(1/Ndata)*abs(fft(xwin,res)).^2;
%pwrSpecXwinApprox=(1/Ndata)*abs(fft(ewin,res)).^2 .* abs(fft(h,res)).^2;
pwrSpecXwinApprox=abs( fft(xwinApprox,res)).^2 / length(xwinApprox);


% generate plots for windowed input convolved with filter and windowed
% output to display approximation errors
figure(1)
clf(figure(1))
title('Windowed output vs Conv. of windowed input')
subplot(4,1,1);
plot(h);
title('filter')

subplot(4,1,2);
plot(x);
hold on
plot([(ceil(N/2)+1):N],xwin,'r')
title('filter output, red part is windowed output for PSD analysis')

subplot(4,1,3);
plot(e);
hold on
plot([(ceil(N/2)+1):N-L+1],ewin,'r')
title('filter input, red part is windowed input for PSD analysis')


subplot(4,1,4)
plot(xwin)
hold on
plot(xwinApprox,'r')
title('windowed filter output (blue), filtered windowed input (red)')



figure(2)
clf(figure(2))
title('PSD of Xwin (blue), PSD of Ewin times abs squared of DTFT of filter (red)');
plot([0:res/2-1]*(1/res),10*log10(pwrSpecXwin(1:res/2)))
hold on
plot([0:res/2-1]*(1/res),10*log10(pwrSpecXwinApprox(1:res/2)),'r')
xlabel('\times 2\pi/T')
ylabel('Magnitude in dB (10*log10(|PSD|))')









