Multiple Power Spectral Density spectra on one plot
13 views (last 30 days)
Show older comments
Hello,
I have calculated the PSD spectrum of a signal and displayed it from 0-500Hz using the following code:
periodogram(M3A,[],'onesided',512,1000)
I would like to display the spectrum up to 2000Hz, and then plot another spectrum on the same graph for comparison.
Please could someone point me in the right direction?
Best Regards
Andy
0 Comments
Accepted Answer
Wayne King
on 15 Dec 2011
Hi Andy, sorry now I understand. You can do this:
% I'll just create an example
Fs = 44100;
t = 0:1/Fs:1-1/Fs;
M3AnO = cos(2*pi*1000*t)+randn(size(t));
M3A = cos(2*pi*500*t)+filter(1,-0.9,rand(size(t)));
[Pxx1,Fxx] = periodogram(M3AnO,[],length(M3AnO),Fs);
[Pxx2,Fxx] = periodogram(M3A,[],length(M3A),Fs);
plot(Fxx,10*log10(Pxx1)); hold on;
plot(Fxx,10*log10(Pxx2),'r'); hold on;
set(gca,'xlim',[0 2000]);
legend('M3AnO','M3A','Location','NorthEast');
xlabel('Hz'); ylabel('dB/Hz');
More Answers (1)
Wayne King
on 15 Dec 2011
Hi Andy, from the code you have shown above, your sampling frequency is 1 kHz. That means your spectrum is periodic with a period equal to 1 kHz. It then sounds like you want to compare it against the power spectral density of another process with a higher sampling rate?
You don't gain anything by extending your process spectrum to 2 kHz because it's periodic. You're better off (if the other process has a higher sampling rate) limiting the part of the spectrum you consider in the process with the higher rate.
Also, the sampling interval and the length of the input time series are part of the scaling factors for the periodogram.
And both the length and sampling interval affect frequency resolution (spacing of the DFT bins).
2 Comments
Ken Saw
on 24 Mar 2014
load sunspot.dat
year=sunspot(:,1);
relNums=sunspot(:,2);
plot(year,relNums)
title('Sunspot Data')
Y = fft(relNums);
Y(1)=[];
n=length(Y);
power = abs(Y(1:floor(n/2))).^2;
nyquist = 1/2;
freq = (1:n/2)/(n/2)*nyquist;
plot(freq,power)
xlabel('cycles/year')
title('Periodogram')
period=1./freq;
plot(period,power);
axis([0 40 0 2e+7]);
ylabel('Power');
xlabel('Period (Years/Cycle)');
hold on;
index=find(power==max(power));
mainPeriodStr=num2str(period(index));
plot(period(index),power(index),'r.', 'MarkerSize',25);
text(period(index)+2,power(index),['Period = ',mainPeriodStr]);
hold off;
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!