Power spectrum of time series

I have the following code with time series of a laser signal output. I need to plot its power spectum or its distribution with the wavelengths, any help please?

function DUMB

clc

tspan=linspace(0,1e-6,10000);

y0=[1e-6 0 0];

[T,Y]= ode45(@rate_equation,tspan,y0);

figure;

plot(T(2000:end)*1e9,Y(2000:end,1),'k');xlabel('Time(ns)');ylabel('E')

xlim([400 420])

dt = T(2)-T(1);

fs=1/dt;

Y1 = fftshift(Y(2000:end,1));

Answer 1

[Mathieu NOE]

hello

it seems that your time data is quite long (~10^5 samples ) , and you are doing just a single buffer fft computation ? you end up with a high resolution spectrum, but this is an spectral average over the entire time vector.

I don't know if your spectrum is supposed to be stable with time - if not, maybe you should do a time / frequency analysis using specgram

I suggest this alternative (here for audio processing, but you can easily adapt it to your application)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dummy signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fs = 5000;
samples = 255001;
dt = 1/Fs;

t = (0:dt:(samples-1)*dt);
omega = 2*pi*(25+20*t);

sensordata = randn(size(t)) + 5*sin(omega.*t);  % signal = linear chirp + noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% FFT parameters
samples = length(sensordata);
NFFT = 512;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% option 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sensor_spectrum, freq] = pwelch(sensordata,w,NOVERLAP,NFFT,Fs);

figure(1),plot(freq,20*log10(sensor_spectrum));
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude (dB)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% option 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

saturation_dB = 80;  % dB range scale (means , the lowest displayed level is 80 dB below the max level)

[sg,fsg,tsg] = specgram(sensordata,NFFT,Fs,w,NOVERLAP);  

% FFT normalisation and conversion amplitude from linear to dB peak
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

% saturation to given dB range
mini_dB = round(max(max(sg_dBpeak))) - saturation_dB;
sg_dBpeak(sg_dBpeak<mini_dB) = mini_dB;

% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);axis('xy');colorbar('vert');
title(['Spectrogram  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');

[Hassan Sultan]

Thank you Mathieu. I'll try.