Power spectrum of time series
10 views (last 30 days)
Show older comments
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));
0 Comments
Answers (1)
Mathieu NOE
on 26 Oct 2020
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)');
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!