PDA

View Full Version : COP Power spectral analysis - little help



dcribeiro22
04-05-2011, 02:24 AM
Hi all,

'Kia Ora'!

Greetings from New Zealand!

We are PhD Students at the University of Otago and we want to calculate total power and median power frequency for COP during standing bipedal and unipedal task.

The total power calculation was based on Soeda et al paper (Exp Brain REs, 2003, 148:226-271).

We are having problems to 'cross-check' our results, specially, regarding total power value. The magnitude for total power is 0.41184.
This refers to a bipedal task, with 30 seconds of duration, fs = 400 Hz, COP measured in meters. The data has already been filtered (low pass, butterworth, 4th order, cutt-off frequency 10 Hz).

Currently, we are working with the code described below.
Would anyone be available to check if this code makes sense? Or maybe you may have a better solution?

We appreciate the attention and help.

Kind regards,

Daniel and Ram


------------------------------------
Octave/Matlab code:

close all
clear all
data = load('bipedal_offset.txt');
fs = 400; % Sample frequency (Hz)
m = length(data(:,2)); % Window length
n = 2^(nextpow2(m)); % Transform length
y = fft(data(:,2),n); % DFT
f = (0:n-1)*(fs/n); % Frequency range
p = y.*conj(y)/n; % Power of the DFT

[nl,nc] = size(data);

plot(f,p)
xlabel('Frequency (Hz)')
ylabel('Power')
title('{\bf Periodogram}')

total_power = sum(sqrt(p).^2)

bogert
04-06-2011, 10:03 AM
Hi Daniel,

According to Parseval's theorem (http://en.wikipedia.org/wiki/Parseval%27s_theorem), you do not need a Fourier transform to calculate the total power. You can simply integrate the square of the signal itself. You may want to first subtract the DC component, because the average of the COP is not of interest to you, only the variations. The Matlab code would be:


total_power = sum((data-mean(data)).^2);

Summing the squared signal after subtracting the average, that is exactly like calculating a standard deviation. So you can also compute this with:



total_power = (n-1)*std(data)^2;


Both should give exactly the same result as what you did with the FFT. One thing that seems wrong with your code is that you compute the total squared magnitude of the entire spectrum, including the DC component at p(1). If you are not interested in the DC component (the average COP location), you may need to do:



total_power = sum(sqrt(p(2:end)).^2);


If that gives the same answer as my code, I would say you can trust either method.

Because of the close relationship between toral power and SD, why not simply report the SD? I suppose it depends on whether you are a statistician or an electrical engineer.

If you want to get a variable that is independent of the length (n) of the signal, you may need to use mean power instead of total power (mean_power = total_power/n). The SD already is normalized that way.

Ton van den Bogert

bogert
04-06-2011, 10:10 AM
Oh, and by the way you were doing the square root (sqrt), and then you squared the result (.^2). Although not wrong, it is unnecessary. You can simply do:



total_power = sum(p(2:end));


Ton

bogert
04-06-2011, 10:20 AM
I tested this on a random signal and indeed all three methods give the same answer:



>> n = 100;
>> data = rand(n,1);
>> sum((data-mean(data)).^2)

ans =

8.6627

>> (n-1)*std(data)^2

ans =

8.6627

>> y = fft(data);
>> p = y.*conj(y)/n;
>> sum(p(2:end))

ans =

8.6627

>>

dcribeiro22
04-07-2011, 08:17 AM
Dear Ton,

Thank you for your comments. Indeed we have removed subtract the DC component. We forgot to inform it...

We really appreciate your help!!!!

Now we feel confortable to continue with data analysis.

All the best,

Daniel and Ram

safaee
04-20-2011, 12:15 AM
Hi
There is a GUI in MATLAB (SPTOOL) which is useful for analysing, filtering and making spectrul analysis too with diffrent methods.
Good Luck

dcribeiro22
04-22-2011, 07:45 PM
Hi
There is a GUI in MATLAB (SPTOOL) which is useful for analysing, filtering and making spectrul analysis too with diffrent methods.
Good Luck


Dear Zahra,

Thank you for your message!
We will check that too.

Kind regards