% Test FFT of random noise
% FFTofNoise.m  L Braile 3/18/06
n= [64 128 256 512 1024];   % Number of points in noise signal
dt = 0.01;    % sample interval (fnyq = 50 Hz)
SS = zeros(5,2048);
for i = 1:5
    npts = n(i);
    s = randn(1,npts); % calculate random noise
    nf = 2048;         % set length of fft
    S = fft(s,nf); % calculate FFT of noise, the added zeros (nf = 1024)
% do not change the spectrum (except by providing finer sampling, df)
% because the zeros don't contribute to the sum of the area in the 
% Fourier integral
    fnyq = 1/(2*dt);  % Nyquist frequency
    df = fnyq/(nf/2); % calculate frequency sample interval
    f = [0:df:fnyq];  % calculate frequency variable, will be (nf/2)+1 long
    SS(i,:) = S.*conj(S)/npts;  
    SS(i,:) = sqrt(SS(i,:)); % calculate amplitude spectrum
    mean(SS(i,:))
    std(SS(i,:)) 
end
t = [0:dt:n(5)*dt - dt];
plot(t,s,'r-')
set(gca,'fontsize',16,'linewidth',2)
xlabel('Time (s)','fontsize',16)
ylabel('Amplitude','fontsize',16)
title('Random Noise Time Series','fontsize',16)
figure
plot(f,SS(:,1:(nf/2)+1),'-')
set(gca,'fontsize',16,'linewidth',2)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Amplitude','fontsize',16)
title('Amplitude Spectrum of Noise','fontsize',16)
figure
plot(f,SS(:,1:(nf/2)+1),'-','linewidth',2)
set(gca,'fontsize',16,'linewidth',2)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Amplitude','fontsize',16)
title('Amplitude Spectrum of Noise','fontsize',16)
axis([0 10 0 3]);
figure
semilogy(f,SS(:,1:(nf/2)+1),'-','linewidth',2)
set(gca,'fontsize',16,'linewidth',2)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Amplitude','fontsize',16)
title('Amplitude Spectrum of Noise','fontsize',16)
axis([0 10 0.01 10]);
