% This routine performs a FFT of normally distributed random noise with a sine wave
%
clear all;
t=[0:400]';
r=randn(size(t));

%add a sine wave with period of 25 -- the amplitude is the multiplier before the sine term
r=r+(.3*sin(2*pi*t/25));

%subtract the mean
r=detrend(r);
figure(1);plot(t,r); 
%%

%calculate FFT and Power after padding
npt=2^14;
H=fft(r,npt);
P=H.*conj(H);
%calculate frequency array
Nyquist=0.5/(t(2)-t(1));
f=linspace(0,Nyquist,npt/2);
%normalize to unit mean power
P=P/mean(P(1:npt/2));
%plot
fmax=.3;
num=round(npt/2*fmax/Nyquist);
figure(2);
plot(f(1:num),P(1:num))
xlabel('frequency')
ylabel('spectral power')
mark1xx

%% simple fft

[ P,f ] = dbfft( t,r,.3,2^14);
figure(3); plot(f,P);
xlabel('frequency')
ylabel('spectral power')
mark1xx

%% noise simulation



%misc steps before iteration
steps=1000; %number of iteration steps
Power=zeros(num,steps); %creates an empty array for the spectral power
w=0; % a counter

for i=1:steps
	nr=randn([(length(t)) 1]); % generates a string of random numbers
	w=w+1; %advances the counter
	%calculate FFT and Power after padding
	Hn=fft(detrend(nr),npt)'; % puts the result into a column vector
	Pn=Hn.*conj(Hn); % gets rid of teh imaginary part
	%normalize to unit mean power
	Pn=Pn/mean(Pn(1:npt/2));
	Power(:,w)=Pn(1:num)'; % loads the Power column into the larger array
end
PQ=mean(Power'); % calculates the mean of each column of the transpose of Power
Pmax=max(Power'); % maximum power at each frequency
PR=std(Power');  % calculates the standard deviation of each column of the transpose of Power
 
P2S=3*PQ; % adds 3 standard deviations to the mean

figure(4); plot(f,P,f,P2S,f,Pmax);
legend('data','95% cl','max noise');
xlabel('frequency')
ylabel('spectral power')
mark1xx

%% bandpass

[ rnew ]=sqbandpass( t,r,27,23 );
figure(5); plot(t,r,t,rnew);
legend('data','bandpass')

