Thanks a lot for your response, your code work fine for me.
Now I am getting slope m = -24.1276 and c = 42.9478 of linear fit (y= mx +c) . Where as when I am analyzingthe same data by originlab software I am getting m = 2.41134 and c = 8.0639.
Probably I need to multiply the pdsx with a factor of 100 ??
Could you please tell me what is the explanation of multiplying psdx with factor of 100 ??
"Matthew Crema" wrote in message <email@example.com>... > "Rajan" wrote in message <firstname.lastname@example.org>... > > Hello, > > > > Please help me perform Power Spectral analysis of electrophysiological signal data recorded from an ion channel. I want to take fft of the Signal and plot it with frequency (log-log plot). and then I want linear fitting of the Power spectra. So that I can Calculate Slope or Noise Properties of PDS. > > > > The Sampling frequency is 10 kHz > > Length of signal/Number of data point is 4096 > > Here is the code I am using > > > > Fs = 10000; > > N = length(x); > > xdft = fft(x); > > xdft = xdft(1:N/2+1); > > psdx = (1/(Fs*N)).*abs(xdft).^2; > > psdx = 2*psdx; > > freq = (0:Fs/length(x):Fs/2)'; > > psdx = 10*log10(psdx); > > freq = (log10(freq)); > > plot(freq,psdx); grid on; > > title('Periodogram Using FFT'); > > xlabel('Frequency (10 kHz)'); ylabel('Power (dB/Hz)'); > > > > p = polyfit(freq,psdx,1) > > f = polyval(p,freq); > > hold on > > plot(freq,f,'r') > > > > The Problem is that I am getting completly wrong Slope of my power spectra (I am comparing my results with the results of originlab software analysis.). > > > > Please let me know what is wrong in my code. - > > How far are you off? Nothing jumps out as wrong, except there may be some confusion about units. I believe that your y-axis label is incorrect as what you are plotting has units of power (dB) and not power spectral density (to get dB/Hz you must divide by sample frequency). > > This may provide a clue into why your slope is not what you expect. What units would you like to report the slope in? eg. dB/decade dB/octave dB/Hz? > > If you want dB/decade I think you want to change to > p = polyfit(10*freq, pdsx, 1) > > A side note: if you have the signal processing toolbox, much of your code can be reduced to one line using PERIODOGRAM, reducing the chances that there is an error that I am missing > > Try this: > N = length(x); > [psdx, freq] = periodogram(x, rectwin(N), N, Fs, 'onesided', 'ms'); > freq = log10(freq); > psdx = 10*log10(psdx); > plot(freq,psdx, 'k-'); grid on; > title('Periodogram Using FFT'); > xlabel('Frequency (10 kHz)'); ylabel('Power (dB)'); > > % Ignore first freq point which is -Inf > p = polyfit(freq(2:end),psdx(2:end),1) > f = polyval(p,freq); > hold on > plot(freq,f,'r') > fprintf('Slope = %g dB/dB\n', p(1))