
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Code for DSP Lecture of April 2, 2005
% 
%  Different parts are separated by "pause"
%   You need to press a key to proceed
%   to following part
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%% Power Spectrum Estimation Using Periodogram Averaging

%%%%% White Gaussian Noise generation
clear all
N= 4096;
x= randn(N,1);

%% Method 1    - Direct Method 
fx= fftshift(fft(x));
Pxx1= abs(fx).^2 /N;
figure(1);  title('Pure White Gaussian Noise');
subplot(4,1,1); plot(x); xlabel('time'); ylabel('Amplitude');
subplot(4,1,2); plot(Pxx1);  xlabel('Frequency'); ylabel('Pxx');

%% method 2     - Nonoverlapped Periodogram Averaging
M= 64;                 % You should try different values of M to see its effect
Pxx2= zeros(M,1);
for k=1:N/M-1 
  s= x(k*M:(k+1)*M-1);  
  fs= fftshift(fft(s));
  Pxx2= Pxx2+abs(fs).^2;
end
Pxx2= Pxx2/N;
subplot(4,1,3); plot(Pxx2);  xlabel('Frequency'); ylabel('Pxx');
   
%% method 3      - Overlapped Periodogram Averaging
M= 64;                 % You should try different values of M to see its effect
Pxx3= zeros(M,1);
for k=1:N-M 
  s= x(k:(k+M-1));  
  fs= fftshift(fft(s));
  Pxx3= Pxx3+abs(fs).^2;
end
Pxx3= Pxx3/N/(N-M)*(N/M);
subplot(4,1,4); plot(Pxx3);  xlabel('Frequency'); ylabel('Pxx');

%% notice the range of variation of power spectra - obvious from the Pxx scale.



%%%%%% Same Example with a Sinusoidal Signal Embedded in Noise    
t= 1:N;
x= x + 0.5*sin(2*pi/13.*t');

%% Method 1    - Direct Method 
fx= fftshift(fft(x));
Pxx1= abs(fx).^2 /N;
figure(2);  title('Pure White Gaussian Noise');
subplot(4,1,1); plot(x); xlabel('time'); ylabel('Amplitude');
subplot(4,1,2); plot(Pxx1);  xlabel('Frequency'); ylabel('Pxx');

%% method 2     - Nonoverlapped Periodogram Averaging
M= 64;                 % You should try different values of M to see its effect
Pxx2= zeros(M,1);
for k=1:N/M-1 
  s= x(k*M:(k+1)*M-1);  
  fs= fftshift(fft(s));
  Pxx2= Pxx2+abs(fs).^2;
end
Pxx2= Pxx2/N;
subplot(4,1,3); plot(Pxx2);  xlabel('Frequency'); ylabel('Pxx');
   
%% method 3      - Overlapped Periodogram Averaging
M= 64;                 % You should try different values of M to see its effect
Pxx3= zeros(M,1);
for k=1:N-M 
  s= x(k:(k+M-1));  
  fs= fftshift(fft(s));
  Pxx3= Pxx3+abs(fs).^2;
end
Pxx3= Pxx3/N/(N-M)*(N/M);
subplot(4,1,4); plot(Pxx3);  xlabel('Frequency'); ylabel('Pxx');

%% Notice the difference in frequency resolution and the difficulty to see
%% the sinusoidal signal in the time domain 

%%%%%%%% Linear Prediction
% doc lpc
L= 1024;
LP= 16;
x = SQUARE(2*pi*(1:L)'/32.)+0.1*randn(L,1);
A= lpc(x, LP);
xc = filter([0 -A(2:end)],1,x);    % Estimated signal   
e = x - xc;                        % Prediction error
figure(3);
subplot(3,1,1); plot(x); xlabel('time'); ylabel('Original');
subplot(3,1,2); plot(xc);xlabel('time'); ylabel('Estimated');
subplot(3,1,3); plot(e);xlabel('time'); ylabel('Error');

%% try to use different functions and different LP to see the different
