Power spectrum of time series
18 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Mathieu NOE
le 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)');
Voir également
Catégories
En savoir plus sur Spectral Measurements dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!