Prev: memory usage issue
Next: Publish stylefile
From: Shaun Hurley on 19 Jul 2010 12:35
I am looking to a Fast Fourier Transformation with a set of data. All of the tutorials I have seen so far used a function to describe the data, rather than data points. I am also trying to transform the data into the power spectral density if anyone is familiar with that, but for now I’m trying to take care of one thing at a time. I did my best to try and follow an example (http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.html ). I used a graph digitalizer to pick off points from the graph (shown below as ts and y). I’ll tell you right now that I am fairly new to MATLAB and don’t fully understand why they did everything they did in the example. Any advice you can offer would be really appreciated. Here is what I have so far, please let me know if you have any suggestions, thank you!

clc
clear all

ts=[0.135465 1.01421 1.83959 2.06564 3.02799 4.99798 5.92092 6.95292 8.04732 9.00228 9.73201 10.9037 12.084 12.7356 13.882 15.0719 15.9746 16.9778 18.0098 19.0847 19.9016 20.9272 23.0615 23.9903 25.1048 25.8365 26.7887 27.8954 29.0743 29.9381 31.1815 31.7711 32.9439 34.0212 34.9153 35.9418 37.061 37.8498 38.9357 40.0374 40.9677 41.9211 42.8185 43.193 43.8673 46.0427 46.8459 48.0351 49.0544 ]/1000;
y=[-1.09364 -1.04245 2.8032 4.03281 6.07267 4.81402 -0.682513 0.676884 -0.0310899 0.595698 0.149883 0.226855 1.95241 0.564645 -4.19947 -0.64215 3.98844 -0.173616 1.18578 -3.23805 -1.01482 -0.885325 1.28384 -3.08746 0.0512647 -0.0281973 0.0752297 1.72238 3.16009 0.35896 0.147987 0.906094 1.19241 -2.78657 0.195438 0.508114 4.53656 1.39534 -0.935056 -0.256124 -4.33958 -4.00064 -0.390604 1.23134 4.18746 0.207047 -0.186527 3.23995 2.16572];
% plot(ts,y)

fs = 1./ts;
L = length(y);
b = 1:L;
t = 50/1000;
% ts = Ts*(b-1)

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L; % fft transformation of y
f = fs./2.*linspace(0,1,NFFT/2+1);

plot(f,2*abs(Y(1:NFFT/2+1)));
From: Wayne King on 19 Jul 2010 13:12
"Shaun Hurley" <jumpinghurley(a)yahoo.com> wrote in message <i21usa$j1n$1(a)fred.mathworks.com>...
> I am looking to a Fast Fourier Transformation with a set of data. All of the tutorials I have seen so far used a function to describe the data, rather than data points. I am also trying to transform the data into the power spectral density if anyone is familiar with that, but for now I&#8217;m trying to take care of one thing at a time. I did my best to try and follow an example (http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.html ). I used a graph digitalizer to pick off points from the graph (shown below as ts and y). I&#8217;ll tell you right now that I am fairly new to MATLAB and don&#8217;t fully understand why they did everything they did in the example. Any advice you can offer would be really appreciated. Here is what I have so far, please let me know if you have any suggestions, thank you!
>
> clc
> clear all
>
> ts=[0.135465 1.01421 1.83959 2.06564 3.02799 4.99798 5.92092 6.95292 8.04732 9.00228 9.73201 10.9037 12.084 12.7356 13.882 15.0719 15.9746 16.9778 18.0098 19.0847 19.9016 20.9272 23.0615 23.9903 25.1048 25.8365 26.7887 27.8954 29.0743 29.9381 31.1815 31.7711 32.9439 34.0212 34.9153 35.9418 37.061 37.8498 38.9357 40.0374 40.9677 41.9211 42.8185 43.193 43.8673 46.0427 46.8459 48.0351 49.0544 ]/1000;
> y=[-1.09364 -1.04245 2.8032 4.03281 6.07267 4.81402 -0.682513 0.676884 -0.0310899 0.595698 0.149883 0.226855 1.95241 0.564645 -4.19947 -0.64215 3.98844 -0.173616 1.18578 -3.23805 -1.01482 -0.885325 1.28384 -3.08746 0.0512647 -0.0281973 0.0752297 1.72238 3.16009 0.35896 0.147987 0.906094 1.19241 -2.78657 0.195438 0.508114 4.53656 1.39534 -0.935056 -0.256124 -4.33958 -4.00064 -0.390604 1.23134 4.18746 0.207047 -0.186527 3.23995 2.16572];
> % plot(ts,y)
>
> fs = 1./ts;
> L = length(y);
> b = 1:L;
> t = 50/1000;
> % ts = Ts*(b-1)
>
> NFFT = 2^nextpow2(L); % Next power of 2 from length of y
> Y = fft(y,NFFT)/L; % fft transformation of y
> f = fs./2.*linspace(0,1,NFFT/2+1);
>
> plot(f,2*abs(Y(1:NFFT/2+1)));

Hi Shaun, Is ts the vector of time measurements from your data? It appears so from your code example, but if this is accurate your data are not uniformly sampled, i.e. the sampling interval is not the same between measurements.

Wayne
From: Shaun Hurley on 19 Jul 2010 13:22
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.
From: Ravi on 19 Jul 2010 13:37
Hi Shaun,

Examples often have extraneous stuff that may not apply to all situations. Below is the minimum usage required.

Basically MATLAB fft function does a discrete Fourier transform for a vector data. This is what applies to your situation. The basic function call is

FT = fft(y,N);

y=vector data at ny points
N=positive integer N>=ny
FT=vector of N complex Fourier coefficients

Often N is not specified : FT = fft(y) then N=ny by default.

In the MATLAB example, N is the power of 2 just higher than ny (obtained using fnc=nextpow2). This optimizes the algorithm, but you don't have to use it. This will also increase the number of Fourier coeff.

The function call is quite simple. However interpretation of the results is often not simple, depending on what you are trying to do. Knowledge of Fourier transformation is essential. The following points are important.

In traditional Fourier series, a factor of 1/N is applied to each coefficient.
In MATLAB this factor is not applied to coeff but to the overall series.
So if you are using the Fourier series, you should apply a factor of 1/N to MATLAB coeff FT.

Fourier coefficients are complex - so they occur in pairs, except for the central (Nyquist) term. Often you have to extract one-sided coefficients, such as for plotting magnitudes. This depends on whether N is odd or even.

You can use the following code to extract one-sided coefficients.

if (rem(N,2)==0)
% even no. points
L=N/2;
else
% odd no. points
L=(N-1)/2;
end

lb=1;
ub=L+1;

FT_onesided=FT(lb:ub) % one-sided coeff

This odd/even scenario should be considered if you are planning to use Fourier series to obtain a harmonic representation to the data y, but by using only a subset of N terms (say a p-order series). A little more processing is required for this, using the above odd/even logic.

Hope this helps !

Ravi





"Shaun Hurley" <jumpinghurley(a)yahoo.com> wrote in message <i21usa$j1n$1(a)fred.mathworks.com>...
> I am looking to a Fast Fourier Transformation with a set of data. All of the tutorials I have seen so far used a function to describe the data, rather than data points. I am also trying to transform the data into the power spectral density if anyone is familiar with that, but for now I&#8217;m trying to take care of one thing at a time. I did my best to try and follow an example (http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fft.html ). I used a graph digitalizer to pick off points from the graph (shown below as ts and y). I&#8217;ll tell you right now that I am fairly new to MATLAB and don&#8217;t fully understand why they did everything they did in the example. Any advice you can offer would be really appreciated. Here is what I have so far, please let me know if you have any suggestions, thank you!
>
> clc
> clear all
>
> ts=[0.135465 1.01421 1.83959 2.06564 3.02799 4.99798 5.92092 6.95292 8.04732 9.00228 9.73201 10.9037 12.084 12.7356 13.882 15.0719 15.9746 16.9778 18.0098 19.0847 19.9016 20.9272 23.0615 23.9903 25.1048 25.8365 26.7887 27.8954 29.0743 29.9381 31.1815 31.7711 32.9439 34.0212 34.9153 35.9418 37.061 37.8498 38.9357 40.0374 40.9677 41.9211 42.8185 43.193 43.8673 46.0427 46.8459 48.0351 49.0544 ]/1000;
> y=[-1.09364 -1.04245 2.8032 4.03281 6.07267 4.81402 -0.682513 0.676884 -0.0310899 0.595698 0.149883 0.226855 1.95241 0.564645 -4.19947 -0.64215 3.98844 -0.173616 1.18578 -3.23805 -1.01482 -0.885325 1.28384 -3.08746 0.0512647 -0.0281973 0.0752297 1.72238 3.16009 0.35896 0.147987 0.906094 1.19241 -2.78657 0.195438 0.508114 4.53656 1.39534 -0.935056 -0.256124 -4.33958 -4.00064 -0.390604 1.23134 4.18746 0.207047 -0.186527 3.23995 2.16572];
> % plot(ts,y)
>
> fs = 1./ts;
> L = length(y);
> b = 1:L;
> t = 50/1000;
> % ts = Ts*(b-1)
>
> NFFT = 2^nextpow2(L); % Next power of 2 from length of y
> Y = fft(y,NFFT)/L; % fft transformation of y
> f = fs./2.*linspace(0,1,NFFT/2+1);
>
> plot(f,2*abs(Y(1:NFFT/2+1)));
From: Wayne King on 19 Jul 2010 13:52
"Shaun Hurley" <jumpinghurley(a)yahoo.com> wrote in message <i221ju$elk$1(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.

Fs = 1000;
NFFT = 64;
ts = linspace(0,(length(y)*(1/Fs)),length(y));
% 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;
PSDEst(2:end-1)=2*PSDEst(2:end-1);
freq = Fs/2.*linspace(0,1,NFFT/2+1);
figure;
plot(freq,10*log10(PSDEst)); grid on;

Hope that helps,
Wayne
 |  Next  |  Last
Pages: 1 2
Prev: memory usage issue
Next: Publish stylefile