clear
%PSD Test Program Demonstrating Concepts with a White Noise Spectrum
%set a sampling interval (s)
delta=1;

%sampling rate (Hz)
fs=1/delta;

%data length in samples
N=2^15;
t = (0:N-1)*delta;

%synthetic white noise signal, units U
x=(1e-6)*randn(N,1);
figure(1)
plot(t,x)
title ('Time Series')
ylabel('U')
xlabel('t (s)')

%synthetic white noise signal with a 1 Hz sinusoid added for demo purposes
%x=(1e-6)*(randn(N,1)+sin(2*pi*1*(0:N-1)'*delta));

disp('Time Domain Estimates')

%integrated signal energy
p1=sum(x.*x)*delta;
disp(['total signal energy: ',num2str(p1),' (U^2)*s'])

%signal power
p2=p1/(N*delta);
disp(['signal power:  ',num2str(p2),' (U^2)'])

%power spectral density for a white process is equally distributed over the Nyquist interval
p3=p2/(fs/2);
disp(['1-sided signal power per Hz:  ',num2str(p3),' (U^2/Hz)'])
p3db=10*log10(p3);
disp(['1-sided signal power per Hz:  ',num2str(p3db),' (dB rel. 1 U^2/Hz)'])

disp(sprintf('\n'))
disp('Frequency Domain Estimates')

%estimates with pmtm
[Pxx,F]=pmtm(x,4,[],fs);
mPxx=mean(Pxx);
Pxxdb=10*log10(mPxx);

%estimate with pwelch
[Pxxw,Fw]=pwelch(x,N/8,0,[],fs);
mPxxw=mean(Pxxw);
Pxxdbw=10*log10(mPxxw);

df=F(2)-F(1);
dfw=Fw(2)-Fw(1);
p4=sum(Pxx)*df;
p4w=sum(Pxxw)*dfw;

%total energy 
% disp(['signal energy estimated from PSD (PMTM):  ',num2str(p4*N*delta),' (U^2*s)'])
% ==
% disp(['signal energy estimated from PSD (PMTM):  ',num2str(p4*N*delta),' (U^2*s)'])

%total power integrated from the PSD estimate
disp(['integrated PMTM PSD (power):  ',num2str(p4),' (U^2)'])
disp(['integrated Pwelch PSD (power):  ',num2str(p4w),' (U^2)'])

%power spectral density mean level
disp(['PSD mean PMTM PSD level (power/Hz):  ',num2str(mPxx),' (U^2/Hz)'])
disp(['PSD mean Pwelch PSD level (power/Hz):  ',num2str(mPxxw),' (U^2/Hz)'])

disp(['PSD mean 1-sided signal power per Hz (dB; PMTM):  ',num2str(Pxxdb),' (dB rel. 1 U^2/Hz)'])
disp(['PSD mean 1-sided signal power per Hz (dB; Pwelch):  ',num2str(Pxxdbw),' (dB rel. 1 U^2/Hz)'])

figure(2)
semilogx(F,10*log10(Pxx),'g',F,Pxxdb*ones(N/2+1,1),'ro',Fw,10*log10(Pxxw),'b',F,Pxxdbw*ones(N/2+1,1),'k.');
legend('PMTM','Mean Value (PMTM)','Pwelch','Mean Value (Pwelch)','location','best')
title ('1-sided time series PSD')
ylabel ('dB rel. 1 U^2/Hz')
xlabel ('Hz')
axis tight


