Prev: memory usage issue
Next: Publish stylefile
From: Greg Heath on 20 Jul 2010 10:16
On Jul 19, 1:52 pm, "Wayne King" <wmkin...(a)gmail.com> wrote:
> "Shaun Hurley" <jumpinghur...(a)yahoo.com> wrote in message <i221ju$el...(a)fred.mathworks.com>...
> > You're right, ts is the time vector. The samples should be
uniformly spaced out, but the digitizer I used wasn't
perfect. For the sake of making things easier, I would be
perfectly fine with making x = 1 - 49. I changed ts to
ts = linspace(1,49,49), Same problem though.
>
> OK Shaun, let's set aside the issue of the non-uniform
> samples for the moment. I'll just assume your data are
> sampled at 1,000 Hz.

close all, clear all, clc

Fs = 1000; % Sampling rate

> NFFT = 64;
> ts = linspace(0,(length(y)*(1/Fs)),length(y));

% Not quite.

N = length(y)
dts = 1/Fs % Temporal sample spacing
T = N*dts % Temporal period (see below)

ts = linspace(0,T-dts,N);
ts = 0:dts:T-dts;
ts = dts*(0:N-1);

NFFT = N % Sufficient when speed is not necessary

Y = fft(y);

df = Fs/N
df = 1/T

f = linspace(0,Fs-df,N);
f = 0:df:Fs-df;
f = df*(0:N-1);

% The equations for the fft and ifft (use the
% command "doc fft" without the quotes) impose
% the following spectral and temporal periodicities
%
% Y(f+Fs) = Y(f);
% yifft = ifft(Y);
% yifft(t+T) = yifft(t);
%
% Since Fs = N*df and T = N*dts, the corresponding
% discrete index vectors have a periodicity of N.

> % If you have the Signal Processing Toolbox, you
> can do create a PSD estimate for
> % your uniformly sampled data this way:
> H = spectrum.periodogram;
> plot(psd(H,y,'Fs',Fs,'NFFT',NFFT))
> % compare the figure to the "brute" force way
> yDFT = fft(y,NFFT);
> PSDEst = 1/(Fs*length(y))*abs(yDFT(1:33)).^2;

% NDFFT/2 + 1 = 33

> PSDEst(2:end-1)=2*PSDEst(2:end-1);

% To account for discarded negative frequencies)
> freq = Fs/2.*linspace(0,1,NFFT/2+1);
> figure;
> plot(freq,10*log10(PSDEst)); grid on;

This uses a scale factor of (1/Fs) that is not in
the fft documentation power spectrum example. The
reason becomes apparent from the identity (search
"Parceval")

sum(y.^2) = sum(abs(Y).^2/N)

The total power in the waveform y (defined
w.r.t. the continuous time power dissipated
in a 1-ohm resistor)

P = dt*sum(y.^2)
P = sum(abs(Y).^2/(N*Fs))
P = sum(PSD)

Hope this helps.

Greg
First  |  Prev  | 
Pages: 1 2
Prev: memory usage issue
Next: Publish stylefile