% White noise demo for ECEn 370
% Brian Jeffs, 8/25/05
% Respond to pauses by striking any key.

% Generate white noise samples, Hamming window, modified
% periodogram, and averaged periodogram

N = 1e5;
L = 256;
R = L/2;
x = randn(N,1);

w = hamming(L);
v1 = w .* x(1:L);

plot(x(1:L));
title([num2str(L),' samples of white noise x[n]']);
sound(x./max(abs(x)),22000);
pause
KK = input('How many successive windows of data should I display? ');
for k = 0:KK-1
	plot(real(fft(x(1+k*L:L+k*L))));
	title(['Real part of ',num2str(L),' point FFT, X[k]']);
	pause
end

K = fix(N/L);
Xave = zeros(L,1);
for k = 0:K-1
    Xave = Xave + fft(x(1+k*L:L+k*L));
end
Xave = Xave/K;
plot(1:L,real(fft(x(1:L))),1:L,real(Xave));
title(['Real part of ',num2str(L),' point window FFT and averaged FFT for ',num2str(K),' windows']);
legend('first window','averaged FFT');
pause

Pxx = psd(x,L,1,hamming(L),R);
plot(Pxx);
limits = axis;
limits(3)=0;
axis(limits);
title(['Averaged periodogram, L= ',num2str(L),' K = ',num2str(floor(N/(L-R)))]);
pause


% Repeat with LP filtered noise

L = 256;
R = L/2;
w = hamming(L);
v1 = w .* x(1:L);

h = remez(200,[0 .2 .22 .42 .44 1.0],[1 1 .5 .5 0 0]);
[H,W] = freqz(h,1,256);
plot(abs(H));
title('filter frequency response')
pause

x = filter(h,1,x);
x(1:220) = [];  % eliminate filter preload samples.

plot(x(1:L));
title([num2str(L),' samples of filtered x[n]']);
sound(x./max(abs(x)),22000);

pause
plot(real(fft(x(1:L))));
title(['Real part of ',num2str(L),' point FFT']);
pause

K = fix(length(x)/L);
Xave = zeros(L,1);
for k = 0:K-1
    Xave = Xave + fft(x(1+k*L:L+k*L));
end
Xave = Xave/K;
plot(1:L,real(fft(x(1:L))),1:L,real(Xave));
title(['Real part of ',num2str(L),' point window FFT and averaged FFT for ',num2str(K),' windows']);
legend('first window','averaged FFT');
pause

Pxx = psd(x,L,1,hamming(L),R);
plot(Pxx);
limits = axis;
limits(3)=0;
axis(limits);
title(['Averaged periodogram, L= ',num2str(L),' K = ',num2str(floor(N/(L-R)))]);

% Repeat with narrowband BP filtered noise

L = 256;
R = L/2;
w = hamming(L);
v1 = w .* x(1:L);

h = remez(200,[0 .05 .07 .09 .11 1.0],[0 0 1 1 0 0]);
[H,W] = freqz(h,1,256);
plot(abs(H));
title('filter frequency response')
pause

x = filter(h,1,x);
x(1:220) = [];  % eliminate filter preload samples.

plot(x(1:L));
title([num2str(L),' samples of filtered x[n]']);
sound(x./max(abs(x)),22000);

pause
plot(real(fft(x(1:L))));
title(['Real part of ',num2str(L),' point FFT']);
pause

K = fix(length(x)/L);
Xave = zeros(L,1);
for k = 0:K-1
    Xave = Xave + fft(x(1+k*L:L+k*L));
end
Xave = Xave/K;
plot(1:L,real(fft(x(1:L))),1:L,real(Xave));
title(['Real part of ',num2str(L),' point window FFT and averaged FFT for ',num2str(K),' windows']);
legend('first window','averaged FFT');
pause

Pxx = psd(x,L,1,hamming(L),R);
plot(Pxx);
limits = axis;
limits(3)=0;
axis(limits);
title(['Averaged periodogram, L= ',num2str(L),' K = ',num2str(floor(N/(L-R)))]);