"edward kabanyas" <djuky_hmi@yahoo.com> wrote in message <j1crrk$58h$1@newscl01ah.mathworks.com>...
> "Wayne King" <wmkingty@gmail.com> wrote in message <j1b7k8$b7h$1@newscl01ah.mathworks.com>...
> > "edward kabanyas" <djuky_hmi@yahoo.com> wrote in message <j1amkg$clj$1@newscl01ah.mathworks.com>...
> > > Hi Friends,
> > >
> > > I need your help again regarding the power spectrum estimation.
> > >
> > > I have hourly average data along several longitude, for example the data is 24 x 100 matrix, 24 is hour with 1 increment and 100 is longitude step is 1 degree increment. I want to calculate the power spectrum of this data.
> > >
> > > The final result that I want is:
> > > a contour plot of power spectrum with x axis is longitude and yaxis is frequency (1/hour) .
> > >
> > > I have tried to calculate the power spectrum for each longitude (example the first longitude: xxx = data(:,1)) by the following code:
> > >
> > > ts = 1; % time step
> > > t = [0:ts:23]; %hourly
> > > fs = 1/ts;
> > >
> > > m = length(t); % Window length
> > > n = pow2(nextpow2(m)); % Transform length
> > > y = fft(xxx,n); % DFT
> > > f = (0:n1)*(fs/n); % Frequency range
> > > power = y.*conj(y)/n; % Power of the DFT
> > >
> > > After that, I join all "power" form the fisrt longitude until the last one, to make a contour plot in function of frequency. However, it seems the result is not correct. Do you have experience with power spectral density analysis ? Thanks for share !
> > >
> > > Edward
> >
> > Hi Edward, The scaling on the periodogram should not be 1/n where n is the padded length. You should use the length of the data vector.
> >
> > Also, I'm assuming your data is realvalued so you don't need to have the frequency axis go out to 1 cycle/hour. You should just use 1/2 the data and have the frequency axis go from 0 cycles/hour to 1/2 cycle/hour which is your Nyquist frequency.
> >
> > I'll make a matrix the size of yours with noise and then add two sine waves with frequencies of 1/8 cycles/hour and 1/4 cycles/hour at the 10th and 15th longitude.
> >
> > n = 0:23;
> > freq = 0:1/128:1/2;
> > x = 0.05*randn(24,100);
> > x(:,10) = cos(pi/4*n)'+0.01*randn(24,1);
> > x(:,15) = cos(pi/2*n)'+0.01*randn(24,1);
> > xdft = fft(x,128);
> > psd_est = 1/length(x)*abs(xdft(1:65,:)).^2;
> > psd_est(2:end1,:) = 2*psd_est(2:end1,:);
> > contour(1:100,freq,10*log10(psd_est));
> > xlabel('Longitude'); ylabel('Cycles/Hour');
> >
> > Hope that helps,
> > Wayne
>
> Dear Wayne,
>
> Thanks very much for your help. It seems work for my case. However, let me confirm a few things. If we refer to the matlab website (http://www.mathworks.com/help/toolbox/signal/dspdata.psd.html)., they calculate the power spectral density as (only an example):
>
> 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;
>
> % Create a singlesided spectrum
> Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2),'Fs',Fs);
> plot(Hpsd);
>
> From your example, you use Nfft = 128, if we followed the above example, your nfft must be 32 ( nfft = 2^nextpow2(24)). And then should not we divide the power by Fs as this example (Pxx = abs(fft(x,nfft)).^2/length(x)/Fs).
>
> Finally, should not we use fft(x,128, dim) for array data or is it enough fft(x,128) ?
>
> Again, thanks for nice help.
>
> Edward
Hi Edward,
1.) Your sampling frequency is 1 in this case, 1 sample/hour. That is why no Fs appears in the example I wrote. If you had a sampling frequency other than unity, it would appear here:
psd_est = 1/(Fs*length(x))*abs(xdft(1:65,:)).^2;
2.) For matrices, fft() works by default on the columns of the matrix, so you do not have to specify dim in this case. So for a matrix
fft(x) is equivalent to fft(x,[],1)
3.) I think people use the zero pad argument too often for no reason. You don't have to zero pad in many cases. However, in YOUR case you only have 24 data points in your time series. So I think in this case it does help to pad a bit. I just picked a couple values and tried to see which one looked the best without going too far. If you have a priori knowledge that a particular frequency, or frequencies, might be present in the data, then by all means you should pick the zero pad so that the frequencies you expect fall as close as possible to multiples of Fs/N, the spacing of the DFT bins.
Wayne
