Matlab Power Spectrum Plot

2
votes

I want to make a plot of the power spectrum of a particular .wav sound file, over the frequency range from -2000 to +2000 Hz.

Attempt:

This is my code so far:

[s, Fs] = wavread('chord.wav');
Hs=spectrum.periodogram;
psd(Hs,s,'Fs',Fs)

I tried using periodogram algorithm. But the resulting plot is ranged from 0 to 20,000 Hz. So, how can I modify the code so that it would be plotted over -2000 to +2000 Hz instead?

Any help would be greatly appreciated.

2
matlabsignal-processingfft
can a frequency be negative? If you just want to modify the axis range take a look here – shamalaia
you can either modify the axis range as A_C suggested, or use a bandpass filter. – GameOfThrows
What is spectrum.periodogram? What is psd()? I don't think those are standard MATLAB functions/classes? – Tom
@A_C Yes, I just want to modify the axis range. I don't know how to use the syntax in your ink. I added Hpsd = dspdata.psd(s,[-2000:2000]) but I am getting errors. What do I need to do? – Merin

2 Answers

0
votes

I modified one of examples in this support page

    Fs = 32e3;   
    t = 0:1/Fs:2.96;
    x = cos(2*pi*t*1.24e3)+ cos(2*pi*t*10e3)+ randn(size(t));
    nfft = 2^nextpow2(length(x));
    Pxx = abs(fft(x,nfft)).^2/length(x)/Fs;
    Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2),'Fs',Fs);  
    plot(Hpsd)

enter image description here

figure
plot(Hpsd)
axis([9.8 10.2 -90 0])  %change axis range

enter image description here

0
votes

I wrote this code for taking FFT of any time history vector in which T is the time step vector and S is the related displacement or any given quantity you hope to plot their spectrum. ps:STfile is my model displacement time history for all degrees of freedom which I saved it as a .mat file 

    clc
clear

load STfile

% is your starting vector of time

data = S(:,12); % vector of data you want to resample 

N = length(T);

for ii=1:(N-1)
    DT(ii) = T(ii+1)-T(ii);
end

DTS = mean(DT);


data_TS = timeseries(data,T); % define the data as a timeseries 

new_T = T(1):DTS:T(end); % new vector of time with fixed dt=1/fs

data_res = resample(data_TS,new_T); % data resampled at constant fs

plot(data_res)

y=getdatasamples(data_res,1:N);

Fs=1./DTS;
NFFT = 2^nextpow2(N);
Y = fft(y,NFFT)/N; % the division by N is to scale the amplitude

freq = Fs/2*linspace(0,1,NFFT/2+1); % vector of frequencies for the plot
figure()
plot(freq,2*abs(Y(1:NFFT/2+1)),'r') % multiplicated by 2 to recover the energy related
% to the negative frequencies since in this way only half spectrum is plotted
xlabel('Frequency(rad/sec)')
xlim([0 0.05])

example result for this code(time history of heave velocity)