psmej

 Power spectrum code from Mike Jones (6 March 2014)

0001 % Power spectrum code from Mike Jones (6 March 2014)
0002 
0003 function [x yp] = psmej( d , doplot)
0004 %returns power spectrum in K^2/Hz
0005 % for amplitude spectrum plot sqrt(smooth(yp))
0006 %   Detailed explanation goes here
0007    Fs = 100;    % hardwired for CBASS data but pass this as a parameter if you like
0008    L = length(d);
0009    NFFT = 2^nextpow2(L);
0010    d = d-mean(d);
0011    y = fft(d,NFFT);
0012    y = y.*conj(y)/(L*Fs);
0013    yp = (y(1:NFFT/2+1));
0014    x = Fs/2*linspace(0,1,NFFT/2+1);
0015    if (doplot == 1)
0016      semilogy(x,sqrt(yp));
0017    else if (doplot == 2)
0018      loglog(x,sqrt(yp));
0019    end
0020 
0021 end