Multiple Power Spectral Density spectra on one plot

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

Answer 1

[Wayne King]

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');

[Andy]

Thanks Wayne, This is exactly what i needed.

Answer 2

[Wayne King]

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).

[Andy]

Hi Wayne,

The sampling frequency of both signals is 44100. I think i must be reading the demo page wrong!

By changing the 1000 to 44100, it produces the frequency range i require.

I have tried the following code to produce both results on the same graph...

[pxNo, freq] = periodogram(M3AnO);

[px, freq] = periodogram(M3A);

plot(freq, px, freq, pxNo)

However, there doesnt appear to be two plots as i think one may be obscuring the other. The frequency also ranges between 0 and 3.5 as opposed to 0-2000Hz.

When i plot using psd instead of periodogram, both lines are visible, but the frequency axis still does not show Hz.

Thanks for your help.

Regards

Andy

[Ken Saw]

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;