3

I'm using the pwelch method in matlab to compute the power spectra for some wind speed measurements. So, far I have written the following code as an example:

t = 10800; % number of seconds in 3 hours
t = 1:t; % generate time vector
fs = 1; % sampling frequency (seconds)
A = 2; % amplitude
P = 1000; % period (seconds), the time it takes for the signal to repeat itself
f1 = 1/P; % number of cycles per second (i.e. how often the signal repeats itself every second).
y = A*sin(2*pi*f1*t); % signal

fh = figure(1);
set(fh,'color','white','Units', 'Inches', 'Position', [0,0,6,6],...
    'PaperUnits', 'Inches', 'PaperSize', [6,6]);
[pxx, f] = pwelch(y,[],[],[],fs);
loglog(f,10*(pxx),'k','linewidth',1.2);
xlabel('log10(cycles per s)');
ylabel('Spectral Density (dB Hz^{-1})');

I cannot include the plot as I do not have enough reputation points

Does this make sense? I'm struggling with the idea of having noise at the right side of the plot. The signal which was decomposed was a sine wave with no noise, where does this noise come from? Does the fact that the values on the yaxis are negative suggest that those frequencies are negligible? Also, what would be the best way to write the units on the y axis if the wind speed is measured in m/s, can this be converted to something more meaningful for environmental scientists?

Emma Tebbs
  • 1,457
  • 2
  • 17
  • 29

1 Answers1

3

Your results are fine. dB can be confusing.

A linear plot will get a good view,

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
y = sin(2 * pi * 50 * t);     % 50Hz signal

An fft approach,

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
subplot(1,2,1);
plot(f,2*abs(Y(1:NFFT/2+1))) 
xlabel('Frequency (Hz)') 
ylabel('|Y(f)|')

pwelch approach,

subplot(1,2,2);
[pxx, freq] = pwelch(y,[],[],[],Fs);
plot(freq,10*(pxx),'k','linewidth',1.2);
xlabel('Frequency (Hz)'); 
ylabel('Spectral Density (Hz^{-1})');

enter image description here

As you can see they both have peak at 50Hz.

Using loglog for both,

enter image description here

So "noise" is of 1e-6 and exists in fft as well, and can be ignored.

For your second question, I don't think the axis will change it will be frequency again. For Fs you should use the sampling frequency of wind speed, like if you have 10 samples of speed in one second your Fs is 10. Higher frequencies in your graph means more changes in wind speed and lower frequencies represent less changes for the speed.

Rashid
  • 4,326
  • 2
  • 29
  • 54
  • Thanks for the answer. I understand that Hz means cycles per second, so I can use the values on the x axis to calculate the period of a signal. However, what does the per Hz represent on the yaxis? How can this be used? For example, what does a spectral density of 0.8 Hz-1 actually mean? – Emma Tebbs Nov 23 '14 at 09:01
  • @EmmaTebbs, In `y` axis it means the values of `y` are for one value of frequency. Its like a graph speed vs time where speed is (m/s) and time is `s`. – Rashid Nov 23 '14 at 09:09
  • @EmmaTebbs, I think `.8-1 hz` the highest frequency band for you. So It means very sharp changes in wind speed. – Rashid Nov 23 '14 at 09:23
  • Great answer. Thanks. I would upvote but don't have enough rep points yet. – Emma Tebbs Nov 23 '14 at 09:36
  • @EmmaTebbs, Any time. Good luck. – Rashid Nov 23 '14 at 10:50
  • Dear @Rashid, I have some question on the use of pwelch on a signal, should I ask a question? Could you help. It concerns the multiple application of pwelch on windows of a given signal. That is, I don't want to apply pwelch on the original signal, but on parts of it, in order to take a mean vector and a covariance matrix of the output. Is this reasonable? – nullgeppetto Feb 13 '16 at 20:26
  • @nullgeppetto , Hi, I would help as much as I can, I can't think of any problems of applying `pwelch` on parts of a signal. You can ask a question or mail me ghorbani.rashid@gmail.com, Good luck, – Rashid Feb 14 '16 at 06:10
  • 1
    @Rashid, this is really much appreciated! Thanks a lot! I'll try to finish something urgent first and then I will ask a public question here (it might help someone in the future) :) – nullgeppetto Feb 14 '16 at 13:49
  • @Rashid, hi again! I have put some thought on my problem and what I really need is to apply pwelch on the whole signal. However, except for the average results that pwelch returns, I would also like to measure their variances. If I am not mistaken, pwelch does an internal averaging of the partial results. So, I would also like to compute the variances as well. Do you think it's feasible? Thanks once again! – nullgeppetto Feb 19 '16 at 10:43
  • @nullgeppetto, do you want the variance on the frequency domain? – Rashid Feb 20 '16 at 09:46
  • @Rashid exactly! Since pwelch computes a mean over the windows and returns a mean vectors, I would to take the variance vector as well. – nullgeppetto Feb 20 '16 at 09:50
  • @nullgeppetto, I took a look into `pwelch`, `computeperiodogram` and `computepsd` but I couldn't find the code where it combines the PSDs. It doesn't seem so hard, look into the algorithms and put a `var()` in appropriate place. – Rashid Feb 20 '16 at 10:59
  • Thanks a lot @Rashid! – nullgeppetto Feb 20 '16 at 11:06
  • @Rashid, just a silly question. I have found some implementation of pwelch in the Web (looking for pwelch.m), but where should I look for sure? Can I see/edit the code Matlab executes when I call pwelch(...)? – nullgeppetto Feb 20 '16 at 11:16
  • @nullgeppetto, yes of course, it is what I was trying to say in my previous comment. type `open pwelch` or right click on pwelch in your code and click open. You can edit any matlab code and save it (it is better to do ot with another name) and have your own function. Good Luck. – Rashid Feb 20 '16 at 13:14
  • Hi @Rashid, I just made a post about pwelch with some questions. If you have some time it would be nice taking a look! Thanks in advance! Here is the question: http://stackoverflow.com/questions/35528353/using-pwelch-to-a-set-of-signals-some-questions-matlab – nullgeppetto Feb 20 '16 at 20:05