Power spectrum of time series

8 次查看(过去 30 天)
Hassan Sultan
Hassan Sultan 2020-10-26
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));

回答(1 个)

Mathieu NOE
Mathieu NOE 2020-10-26
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)');

类别

Help CenterFile Exchange 中查找有关 Spectral Measurements 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by