Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
power spectrum for multi array

Subject: power spectrum for multi array

From: edward kabanyas

Date: 3 Aug, 2011 05:34:08

Message: 1 of 4

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 y-axis 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:n-1)*(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

Subject: power spectrum for multi array

From: Wayne King

Date: 3 Aug, 2011 10:24:08

Message: 2 of 4

"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 y-axis 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:n-1)*(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 real-valued 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 10-th and 15-th 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:end-1,:) = 2*psd_est(2:end-1,:);
contour(1:100,freq,10*log10(psd_est));
xlabel('Longitude'); ylabel('Cycles/Hour');

Hope that helps,
Wayne

Subject: power spectrum for multi array

From: edward kabanyas

Date: 4 Aug, 2011 01:15:32

Message: 3 of 4

"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 y-axis 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:n-1)*(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 real-valued 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 10-th and 15-th 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:end-1,:) = 2*psd_est(2:end-1,:);
> 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 single-sided 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

Subject: power spectrum for multi array

From: Wayne King

Date: 4 Aug, 2011 09:40:30

Message: 4 of 4

"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 y-axis 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:n-1)*(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 real-valued 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 10-th and 15-th 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:end-1,:) = 2*psd_est(2:end-1,:);
> > 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 single-sided 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

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us