How to calculate power spectrum of an impulse response?
Srdjan
First, I would like to apologize for the lack of my knowledge about Matlab, but I am a biologist and I never had a chance to learn more about Matlab, signal processing and some other stuff that I need now. Currently I doing some electrophysiological experiments and I need to analyze the data. I am measuring responses of photoreceptor cells to brief flashes of light. The cells respond with a short change in voltage potential, which we call an impulse response. Here is an example:
http://img39.imageshack.us/img39/9971/impulseresponse.jpg
If it matters, sampling rate for this response was 10000 points per sec.
From this response I need to calculate power spectra, but I don't know how to do it. Apparently I need to Fourier transform this and square the resulting amplitude spectrum. Here is an example of a power spectrum that was calculated from data similar to mine.
http://img683.imageshack.us/img683/5659/powerspectra.jpg
The graph above has not one, but five responses at five different temperatures. I just need one curve, not five. If its not complicated, can someone explain me how to do this in Matlab. I was playing with the FFT function, but whatever I got was nowhere near the responses shown in the picture above. I have a basic version of Matlab with no additional toolboxes.
Thank you very much for your help!
Godzilla
When you ask for help, it is beneficial to show what you have already done.
Srdjan
I have an impulse response stored in a vector, lets call it "IR", 2501 data points long.
To get the power spectrum I did this:
PS = abs(fft(IR));
Then I normalize the amplitude:
NPS = PS./max(PS);
Second, I tried to calculate the frequency:
Sampling rate was 10000 points/sec, so I defined a variable T=1/10000
Then I did this:
N = 2501;
n = (0:N-1);
fn = n./(N.*T);
At the end I ploted this on a log scale:
semilogx(fn, NPS);
Here is the outcome:
http://img819.imageshack.us/img819/1880/npsg.jpg
I don't know if this is correct. Also I don't understand why this curve starts somewhere around 0.72, because when I look at the NPS vector values go from 1 and then downwards. Here is how the graph looks when I just plot NPS; semilogx(NPS):
http://img46.imageshack.us/img46/2748/nps2.jpg
Thanks for the help!
TideMan
The main problem is that you have asked for N data in the spectrum,
but you only want the data from 2 to fix(N/2). The first component
has zero frequency (DC component), so when you try plotting it on a
log scale, it can't - try this: log10(0). The data from fix(N/2)+1 to
N are essentially a repeat of the first set of data, so you can
discard them.
So:
N2=fix(N/2);
NPS=PS(2:N2)/max(PS(2:n2));
n = (1:N2-1);
fn = n./(N.*T);
semilogx(fn, NPS);
Greg Heath
> Ok, here is where I am right now:
>
> I have an impulse response stored in a vector, lets call it "IR", 2501 data points long.
> To get the power spectrum I did this:
> PS = abs(fft(IR));
N = 2501 % N is odd
Fs = 1e4
f = (Fs/N)*(0:N-1)';
PSD = abs(fft(IR)).^2/(N*Fs); % Power Spectral Density
% Single-sided spectrum for real IR
fss = f(1:(N-1)/2) % for N odd
PSDss = [PSD(1); 2*PSD(2:(N-1)/2)]; % Note the factor of 2
figure
subplot(211)
plot(fss,PSDss)
subplot(212)
semilogx(fss(2:end),PSD(2:end)) % Must omit the fss = 0 point
> Then I normalize the amplitude:
> NPS = PS./max(PS);
>
> Second, I tried to calculate the frequency:
> Sampling rate was 10000 points/sec, so I defined a variable T=1/10000
> Then I did this:
> N = 2501;
> n = (0:N-1);
> fn = n./(N.*T);
>
> At the end I ploted this on a log scale:
> semilogx(fn, NPS);
>
> Here is the outcome:http://img819.imageshack.us/img819/1880/npsg.jpg
>
> I don't know if this is correct. Also I don't understand why this curve starts
> somewhere around 0.72,
NP(1) must be excluded from the semilogx plot because f(1) = 0.
Therefore the first point is NPS(2) = 0.72.
However, since you tried to plot all 2501 points, you should have
conjugate symmetry about the midpoint of f((N-1)/2) and f((N+1)/2 ) .
In particular, you should have NPS(N) = NPS(2) = 0.72.
Wha hoppen?
> because when I look at the NPS vector values go from 1 and then
> downwards. Here is how the graph looks when I just plot NPS; semilogx
>(NPS):http://img46.imageshack.us/img46/2748/nps2.jpg
The x-axis is just log(index) = log(1:N). Therefore all points are
included.
Also note that NPS(N) = 0.72 , as it should.
Hope this helps.
Greg
Srdjan
Thank you very much for your help. With an addition of squaring to TideMan's code (something I forgot to do in the first place), both of your approaches give the same outcome. Now this makes me confident I have the right calculation to proceed with the analysis, though I still need to learn more about the Fourier Transformation and Power Spectrum.
Thank you very much again!
S.