Calculate PSD and plot smooth curve

Hey all,

I've new to matlab, and have been wrestling with processing and plotting PSD data. I have raw data (time and magnitude in g), as well as already processed data by different software (in freq and g^2/Hz). I'd like to compare PSDs from this software versus whatever I can cook up in matlab.

I need to take the raw data (in csv files, two columns, time and amplitude of g) and run a PSD on it. I'd then like to plot a smoothed version of the PSD to make a profile curve that represents most of the data (i don't want to plot the results of the PSD, too much info).

I've been trying to use different functions in matlab some of which seem obsolete or deprecated, so any help is appreciated, thanks!

JT

Answer 1

[Wayne King]

Hi JT, you can use spectrum.periodogram and spectrum.welch

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+randn(size(t));
plot(psd(spectrum.periodogram,x,'Fs',Fs,'NFFT',length(x)));
figure;
plot(psd(spectrum.welch,x,'Fs',Fs,'NFFT',length(x)));

Answer 2

[Jimes Tooper]

Hey Wayne,

Thanks for this! I was wondering, is there a way to get the plots into different units? I'd prefer g^2/Hz to compare to MIL, SAE, ISO standards etc.. Thank you!

JT

Answer 3

[Wayne King]

Hi, yes, you can just return the output from psd() and plot it as you wish.

psdest = psd(spectrum.periodogram,x,'Fs',Fs,'NFFT',length(x));
plot(psdest.Frequencies, psdest.Data); grid on;
ylabel('g^2/Hz');
xlabel('Hz');

Answer 4

[Jimes Tooper]

Great, thank you!

Last question: in all of the examples I see, including yours, it's not clear to me how I would use my data instead. I figure the following: I have my csv file with two columns, time and g. I import using:

[test] = xlsread('data.csv');

This gives me an array of double (A is any number), but I also want to break up the array, so I get:

x = test(:,1); 
y = test(:,2); 
t = 0:length(x);

I make t just to have a time domain starting at zero instead of an arbitrary second in my data, and then plot x vs y just to make sure it looks like my spectra, it does! :)

Then I go for it:

Fs = 1000; 
h = spectrum.welch; 
plot(psd(h,test,'Fs',fs,'NFFT',t);

The plot that pops up is definitely not what I'm looking for and it calculates in under 2 seconds so I'm wary.

Is this how you would go about using your own data?

Answer 5

[Wayne King]

Hi, you do not want to input test in this example, because test is a matrix. psd() expects a vector input, a 1-D signal. What you want to input is either x or y, whichever is the meaningful time series, also, for 'NFFT', you do not want to input a vector, like t =0:length(x). NFFT should just be a number, the number of FFT bins, use length(x), or length(y).

Finally, make sure that fs is your sampling frequency

'Fs',10000 or 'Fs',1 whatever the sampling frequency of your data is.