Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Help with FFT

Subject: Help with FFT

From: Luca Cerone

Date: 6 Sep, 2010 15:10:10

Message: 1 of 16

Dear all,
I hope some of you can lend me a hand to use well fft().
What I need to do is to evaluate frequencies of a real function f(t).
What I have is an array of measurements Y taken on an interval [tstart tend],
and would like to check if there is a "dominant" frequency in the signal,
or if frequencies have some special distribution.

I couldn't really get how to use fft to do this, even reading the documentation.
I've tried to evaluate this for known function (like sin(x) for x from 0 to 10) but didn't get 1/(2*pi) as a result,
so I assume I'm using the function in the wrong way.

Do you have any advice, or minimal example so that I can make things working?
Thanks a lot in advance,
Cheers, -luca
 

Subject: Help with FFT

From: Miroslav Balda

Date: 6 Sep, 2010 15:26:05

Message: 2 of 16

"Luca Cerone" <luca_cerone#_remove_this#@yahoo.it> wrote in message <i6308i$sai$1@fred.mathworks.com>...
...
SNIP
...
Hi Luca,

Try to make
     X = abs(fft(x));
Find the max peak of X and its index:
     [Xmax,I] = max(X);
The frequency correspondig the peak is
     N = length(x);
     fpeak = (I-1)*fs/N;%,
where
fs = sampling frequancy of x(t)

Best regards,

Mira

Subject: Help with FFT

From: Greg Heath

Date: 6 Sep, 2010 15:59:58

Message: 3 of 16

On Sep 6, 11:10 am, "Luca Cerone" <luca_cerone#_remove_th...@yahoo.it>
wrote:
> Dear all,
> I hope some of you can lend me a hand to use well fft().
> What I need to do is to evaluate frequencies of a real function f(t).
> What I have is an array of measurements Y taken on an interval [tstart tend],
> and would like to check if there is a "dominant" frequency in the signal,
> or if frequencies have some special distribution.
>
> I couldn't really get how to use fft to do this, even reading the documentation.
> I've tried to evaluate this for known function (like sin(x) for x from 0 to 10) but didn't get 1/(2*pi) as a result,
> so I assume I'm using the function in the wrong way.
>
> Do you have any advice, or minimal example so that I can make things working?
> Thanks a lot in advance,
> Cheers, -luca

Show your code, results and an explanation of
why you think it is not correct.

Greg

Subject: Help with FFT

From: Luca Cerone

Date: 6 Sep, 2010 16:13:05

Message: 4 of 16

Hi, I think that the result is not correct because I have to rescale frequencies
in the right range.. The solution given by Mira should make the trick,
except that I'd like to see the whole spectrum distribution before to cut the maximum
(want to check if there is more than on peak)

> Show your code, results and an explanation of
> why you think it is not correct.
>
> Greg

Subject: Help with FFT

From: Miroslav Balda

Date: 6 Sep, 2010 16:34:06

Message: 5 of 16

"Luca Cerone" <luca_cerone#_remove_this#@yahoo.it> wrote in message <i633uh$lrm$1@fred.mathworks.com>...
> Hi, I think that the result is not correct because I have to rescale frequencies
> in the right range.. The solution given by Mira should make the trick,
> except that I'd like to see the whole spectrum distribution before to cut the maximum

Hi Luca,
You may see all the results by plotting X:
     plot((0:N/2-1)*fs/N,X(1:N/2);
and judge, whether there are more important peaks. You also may use the function tolpk from FEX:
     www.mathworks.com/matlabcentral/fileexchange/22662
that finds more peaks exceeding a given tolerance.
Best regards,

Mira

Subject: Help with FFT

From: Luca Cerone

Date: 6 Sep, 2010 17:31:05

Message: 6 of 16

Sorry but I still don't understand what I do wrong:
here is the function I've written:

function [freq,p]=positiveFFT(t,y)
% I assume that the number of samples is even
T=max(t)-min(t); %estimate the period;
N=length(t); %evaluate the number of samples;
p = abs(fft(y))/(N/2);%powers
p = p(1:N/2).^2;%powers only in the positive Range;
freq = [0:N/2-1]/T;

It takes as input the vector of times t and the vector of measurements x.
(for the moment I'm assuming the number of samples is even, but I'd like to
get rid of this and create a more generic function).

If I test it using:

t=linspace(10,20,100);
y=sin(2*pi*3*t);
[freq,p]=positiveFFT(t,y);
stem(freq,p)

The results seem wrong to me.
What I'd expect is to have a frequency range from 0 to 5 (but the freq range stops at
4.9) (since T=10 and the number of sample is 100, the frequency of sampling is 10, but I'm only interested in the positive spectrum).
Also I would expect to have a value close to 1 for the p corresponding to freq=3
and all the other values to be almost 0. Instead I get a peak (of about .7) in correspondence of freq=3, but the nearby values are also 0.1, 0.2 almost.

Can you see the mistake in my code? Can you also advice me how to get rid of the
"even number of samples" hypothesis?

Thanks for all your help, it's really appreciated.
Cheers, -luca



> Hi Luca,
> You may see all the results by plotting X:
> plot((0:N/2-1)*fs/N,X(1:N/2);
> and judge, whether there are more important peaks. You also may use the function tolpk from FEX:
> www.mathworks.com/matlabcentral/fileexchange/22662
> that finds more peaks exceeding a given tolerance.
> Best regards,
>
> Mira

Subject: Help with FFT

From: Wayne King

Date: 6 Sep, 2010 18:14:06

Message: 7 of 16

"Luca Cerone" <luca_cerone#_remove_this#@yahoo.it> wrote in message <i638gp$dsj$1@fred.mathworks.com>...
> Sorry but I still don't understand what I do wrong:
> here is the function I've written:
>
> function [freq,p]=positiveFFT(t,y)
> % I assume that the number of samples is even
> T=max(t)-min(t); %estimate the period;
> N=length(t); %evaluate the number of samples;
> p = abs(fft(y))/(N/2);%powers
> p = p(1:N/2).^2;%powers only in the positive Range;
> freq = [0:N/2-1]/T;
>
> It takes as input the vector of times t and the vector of measurements x.
> (for the moment I'm assuming the number of samples is even, but I'd like to
> get rid of this and create a more generic function).
>
> If I test it using:
>
> t=linspace(10,20,100);
> y=sin(2*pi*3*t);
> [freq,p]=positiveFFT(t,y);
> stem(freq,p)
>
> The results seem wrong to me.
> What I'd expect is to have a frequency range from 0 to 5 (but the freq range stops at
> 4.9) (since T=10 and the number of sample is 100, the frequency of sampling is 10, but I'm only interested in the positive spectrum).
> Also I would expect to have a value close to 1 for the p corresponding to freq=3
> and all the other values to be almost 0. Instead I get a peak (of about .7) in correspondence of freq=3, but the nearby values are also 0.1, 0.2 almost.
>
> Can you see the mistake in my code? Can you also advice me how to get rid of the
> "even number of samples" hypothesis?
>
> Thanks for all your help, it's really appreciated.
> Cheers, -luca
>
>
>
> > Hi Luca,
> > You may see all the results by plotting X:
> > plot((0:N/2-1)*fs/N,X(1:N/2);
> > and judge, whether there are more important peaks. You also may use the function tolpk from FEX:
> > www.mathworks.com/matlabcentral/fileexchange/22662
> > that finds more peaks exceeding a given tolerance.
> > Best regards,
> >
> > Mira

Hi Luca,

you need to have one more sample in your frequency vector and you need to take one more index from the DFT of your signal y. The first bin in the DFT is DC (zero frequency). The Nyquist is N/2+1 (N is even in your case). That's why you're not getting to the Nyquist. Your frequency vector should be length 51 and you should take the first 51 indices of the DFT of your signal.

Fs = 10;
N = length(y);
freq = 0:(Fs/N):5;
ydft = abs(fft(y)).^2;
ydft =(1/N)* ydft(1:N/2+1);
ydft(2:end)=2*ydft(2:end);
plot(freq,ydft);


Wayne

Subject: Help with FFT

From: Wayne King

Date: 6 Sep, 2010 19:19:05

Message: 8 of 16

"Wayne King" <wmkingty@gmail.com> wrote in message <i63b1e$nnk$1@fred.mathworks.com>...
> "Luca Cerone" <luca_cerone#_remove_this#@yahoo.it> wrote in message <i638gp$dsj$1@fred.mathworks.com>...
> > Sorry but I still don't understand what I do wrong:
> > here is the function I've written:
> >
> > function [freq,p]=positiveFFT(t,y)
> > % I assume that the number of samples is even
> > T=max(t)-min(t); %estimate the period;
> > N=length(t); %evaluate the number of samples;
> > p = abs(fft(y))/(N/2);%powers
> > p = p(1:N/2).^2;%powers only in the positive Range;
> > freq = [0:N/2-1]/T;
> >
> > It takes as input the vector of times t and the vector of measurements x.
> > (for the moment I'm assuming the number of samples is even, but I'd like to
> > get rid of this and create a more generic function).
> >
> > If I test it using:
> >
> > t=linspace(10,20,100);
> > y=sin(2*pi*3*t);
> > [freq,p]=positiveFFT(t,y);
> > stem(freq,p)
> >
> > The results seem wrong to me.
> > What I'd expect is to have a frequency range from 0 to 5 (but the freq range stops at
> > 4.9) (since T=10 and the number of sample is 100, the frequency of sampling is 10, but I'm only interested in the positive spectrum).
> > Also I would expect to have a value close to 1 for the p corresponding to freq=3
> > and all the other values to be almost 0. Instead I get a peak (of about .7) in correspondence of freq=3, but the nearby values are also 0.1, 0.2 almost.
> >
> > Can you see the mistake in my code? Can you also advice me how to get rid of the
> > "even number of samples" hypothesis?
> >
> > Thanks for all your help, it's really appreciated.
> > Cheers, -luca
> >
> >
> >
> > > Hi Luca,
> > > You may see all the results by plotting X:
> > > plot((0:N/2-1)*fs/N,X(1:N/2);
> > > and judge, whether there are more important peaks. You also may use the function tolpk from FEX:
> > > www.mathworks.com/matlabcentral/fileexchange/22662
> > > that finds more peaks exceeding a given tolerance.
> > > Best regards,
> > >
> > > Mira
>
> Hi Luca,
>
> you need to have one more sample in your frequency vector and you need to take one more index from the DFT of your signal y. The first bin in the DFT is DC (zero frequency). The Nyquist is N/2+1 (N is even in your case). That's why you're not getting to the Nyquist. Your frequency vector should be length 51 and you should take the first 51 indices of the DFT of your signal.
>
> Fs = 10;
> N = length(y);
> freq = 0:(Fs/N):5;
> ydft = abs(fft(y)).^2;
> ydft =(1/N)* ydft(1:N/2+1);
> ydft(2:end)=2*ydft(2:end);
> plot(freq,ydft);
>
>
> Wayne

Sorry Luca, I neglected to look at your scaling question. Look at your signal as a function of time.
plot(t,y);

See how in each period it isn't actually oscillating between [-1,1] as you would expect. Now compare what you did to the following:

Fs = 100;
t = linspace(0,3,300);
y = sin(2*pi*3*t);
plot(t,y);
N = length(y);
ydft = fft(y);
ydft = ydft(1:N/2+1);
freq = 0:(Fs/N):50;
plot(freq,2/N*abs(ydft));

Wayne

Subject: Help with FFT

From: Greg Heath

Date: 7 Sep, 2010 02:19:20

Message: 9 of 16


PLEASE DO NOT TOP-POST BY PLACING YOUR REPLIES
ABOVE THE FORMER POST!

On Sep 6, 1:31 pm, "Luca Cerone" <luca_cerone#_remove_th...@yahoo.it>
wrote:
> Sorry but I still don't understand what I do wrong:
> here is the function I've written:
>
> function [freq,p]=positiveFFT(t,y)
> % I assume that the number of samples is even
> T=max(t)-min(t); %estimate the period;

Incorrect

> N=length(t); %evaluate the number of samples;

dt = t(2)-t(1); % or mean(diff(t))

T = N*dt % period of ifft(fft(y))

> p = abs(fft(y))/(N/2);%powers

Incorrect name and scaling. See below

> p = p(1:N/2).^2;%powers only in the positive Range;
> freq = [0:N/2-1]/T;
>
> It takes as input the vector of times t and the vector of measurements x.
> (for the moment I'm assuming the number of samples is even, but I'd like to
> get rid of this and create a more generic function).
>
> If I test it using:
>
> t=linspace(10,20,100);
> y=sin(2*pi*3*t);
> [freq,p]=positiveFFT(t,y);
> stem(freq,p)
>
> The results seem wrong to me.
> What I'd expect is to have a frequency range from 0 to 5 (but the freq range stops at
> 4.9) (since T=10 and the number of sample is 100, the frequency of sampling is 10, but I'm only interested in the positive spectrum).
> Also I would expect to have a value close to 1 for the p corresponding to freq=3
> and all the other values to be almost 0. Instead I get a peak (of about .7) in correspondence of freq=3, but the nearby values are also 0.1, 0.2 almost.
>
> Can you see the mistake in my code?
>Can you also advice me how to get rid of the

In addition to the other errors, f0 = 3 is not an integral
multiple of df = 1/T = 1/(N*dt) = Fs/N

> "even number of samples" hypothesis?
>
> Thanks for all your help, it's really appreciated.
> Cheers, -luca

t = 0:dt:T-dt;
t = dt*(0:N-1);
t = linspace(0,T-dt,N);

Fs = 1/dt
df = Fs/N % or 1/T

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

% Notice the duality (t,dt,T) <=> (f,df,Fs)

pi = dt*y.^2; % Instantaneous power
P = dt*sum(y.^2) % Total power
Pav = (1/T)*dt*sum(y.^2) % Average power (dt/T) = 1/N

Y = fft(y);
absY = abs(Y);
PSD = absY.^2/(N*Fs) % Power Spectral Density

% Parceval's Theorem

% sum(y.^2) = sum(absY.^2/N)

P = sum(PSD) % Total power
Pav = (1/Fs)*df*sum(PSD) % Average power (df/Fs) = 1/N

If y is real

k = [ 2 : N ] ;
Y(k) = conj( Y(N+2-k) ) ;

Number of unique points

Q = 1 + numel( find( k <= N+2-k) )
Q = ceil( (N+1)/2 )

Bipolar frequency interval

fb = df *[ -ceil((N-1)/2) : floor((N-1)/2)];

if mod(N,2)~=0 % N odd
     fb = df*[ -(N-1)/2 : (N-1)/2];
else % N even
     fb = df*[ -N/2 : N/2-1]; % Note Nyquist < 0
end

Singlesided frequency interval

fss = f(1:Q);
fss = df*[0:Q-1];
fss = df*[0 : ceil((N-1)/2 ];

if mod(N,2)~=0 % N odd
     fss = df*[ 0 : (N-1)/2 ];
else % N even
     fss = df*[ 0 : N/2 ]; % Note Nyquist > 0
end

Singlesided amplitudes and Power Spectra

if mod(N,2)~=0 % N odd

    absYss = [ absY(1) 2*absY(2:Q) ];
    PSDss = [ PSD(1) 2*PSD(2:Q) ];

else % N even

    absYss = [ absY(1) 2*absY(2:Q-1) absY(Q)];
    PSDss = [ PSD(1) 2*PSD(2:Q-1) absY(Q)];
end

Note that only the interior points of PSDss are
proportional to absYss.^2;

Hope this helps.

Greg

Subject: Help with FFT

From: Luca Cerone

Date: 7 Sep, 2010 08:29:04

Message: 10 of 16

Hope I'm not top-posting again (sorry but I don't know what it means..
Does it mean I haven't replied to the last post? if so I'm sorry, didn't know there was
such a convention..)
Thanks Wayne and Greeg, your answer helped me understanding a few things
(Greg thanks for all the "side-measures" to evaluate, some of them I really don't remember
what they're, I'll try to look on Wikipedia and make sense of it, it's a few years I don't use fft and signal processing).

Greg Heath <heath@alumni.brown.edu> wrote in message <e80256a7-de5e-40d9-b922-eec65248c4ed@z28g2000yqh.googlegroups.com>...
>
> PLEASE DO NOT TOP-POST BY PLACING YOUR REPLIES
> ABOVE THE FORMER POST!

Subject: Help with FFT

From: Luca Cerone

Date: 7 Sep, 2010 09:11:05

Message: 11 of 16

Hi, there is a thing I can't really understand:
as far as I've understood, in order for the fft to work,
for a periodic signal I have to stop "slightly" before the
end of the period (that is if T is the period I consider only
points from 0 till T-dt).
If I use points from 0 to T I get results that are slightly inaccurate
(this because fft works for periodic function).
In theory though when I measure a signal, I usually have measurements
for it in a range [0,T] and a priori I don't know if is periodic or not.
In theory might be that the signal is periodic (say of period T/2), so
that I'm taking into account 2 periods for my frequency analysis.
Would such a situation lead to wrong results? If so what are the practical "tricks"
to perform a good frequency analysis?

I'm trying to analyse some biological, and would like to check if there
is a dominant frequency and how this change with respect to certain
parameters.
Though the theory of fourier analysis is in principle pretty clear to me,
if I work with functions that are periodic in the [0, 2*pi] range,
I'm having some difficulties (very stupid, I know I'm a bit rusted...) in
dealing with general functions.

Sorry to be wasting so long time on such an easy issue.
Cheers, -luca

Subject: Help with FFT

From: Greg Heath

Date: 8 Sep, 2010 17:28:24

Message: 12 of 16

On Sep 7, 5:11 am, "Luca Cerone" <luca_cerone#_remove_th...@yahoo.it>
wrote:
> Hi, there is a thing I can't really understand:
> as far as I've understood, in order for the fft to work,
> for a periodic signal I have to stop "slightly" before the
> end of the period (that is if T is the period I consider only
> points from 0 till T-dt).
> If I use points from 0 to T I get results that are slightly
> inaccurate (this because fft works for periodic function).

The fft and ifft summation formulae are consistent
with the periodic constraints

 x[i+N] = x[i]

or, equivalently,

x(t+T) = x(t)

where T is the period.

> In theory though when I measure a signal, I usually have
> measurements for it in a range [0,T]

No.

If T is the period, you usually have measurements in the
range [0,tmax]

>and a priori I don't know if is periodic or not.

Yes you do:

If x(tmax) ~= x(0), then T ~= tmax and the fft/ifft
summation formulae are consistent with the constraint
T = tmax+1/Fs.

In fact, if you manually apply the formulas for fft
and ifft you will find that, regardless of the
periodicity of the original x(t), the reconstructed
xr = idft(dft(x))) will be periodic with period
T = N/Fs.

> In theory might be that the signal is periodic
> (say of period T/2), so that I'm taking into
> account 2 periods for my frequency analysis.
> Would such a situation lead to wrong results?

No.

> If so what are the practical "tricks"
> to perform a good frequency analysis?

1. When I want to analyze frequencies of a complicated
function up to f0 Hertz, I try to use a sampling
frequency Fs >= 4*f0.

2. If abs(x[N]-x[1]) is large, I apply windowing
to eliminate the ficticious high frequencies
caused by the discontinuity


Hope this helps.

Greg

Subject: Help with FFT

From: Luca Cerone

Date: 9 Sep, 2010 11:18:06

Message: 13 of 16

Hi Greg, thanks for your reply, but not everything is clear to me..

> >and a priori I don't know if is periodic or not.
>
> Yes you do:
>
> If x(tmax) ~= x(0), then T ~= tmax and the fft/ifft
> summation formulae are consistent with the constraint
> T = tmax+1/Fs.
Sorry I don't get this point. I agree that if x(0)~=x(tmax)
then T ~=tmax, but this doesn't mean that my signal is periodic.

> 1. When I want to analyze frequencies of a complicated
> function up to f0 Hertz, I try to use a sampling
> frequency Fs >= 4*f0.
So, say I'm interested in frequencies between [0 and 0.1] HZ,
my sampling frequency has to be around 0.025??
(say my signal is the solution of an ODE, do I have to choose t
in a way similar to N=1000 (to say I want 1000 points),
t=[0:1/0.025:N/0.025] ?)

> 2. If abs(x[N]-x[1]) is large, I apply windowing
> to eliminate the ficticious high frequencies
> caused by the discontinuity
How can I apply windowing? Is there any function in Matlab?
I've an idea of what windowing does, but only read of it.

Thanks a lot in advance for your help,
Cheers, -Luca

Subject: Help with FFT

From: Greg Heath

Date: 10 Sep, 2010 12:09:52

Message: 14 of 16

On Sep 9, 7:18 am, "Luca Cerone" <luca_cerone#_remove_th...
@yahoo.it> wrote:
> Hi Greg, thanks for your reply, but not everything is
> clear to me..
>
> > >and a priori I don't know if is periodic or not.
>
> > Yes you do:
>
> > If x(tmax) ~= x(0), then T ~= tmax and the fft/ifft
> > summation formulae are consistent with the constraint
> > T = tmax+1/Fs.
>
> Sorry I don't get this point. I agree that if x(0)~=x(tmax)
> then T ~=tmax, but this doesn't mean that my signal is periodic.

Correction: You know if your original signal is not periodic
with period tmax.

However, the fft doesn't care:

Regardless of whether the original function is periodic
or aperiodic, the fft/ifft summation formulae are
consistent with the constraint T = tmax+1/Fs.

Therefore if you have a function with period T0,
the results will be nicer if you can choose
tmax = m*T0-dt (m = 1,2,...; dt = 1/Fs)

> > 1. When I want to analyze frequencies of a complicated
> > function up to f0 Hertz, I try to use a sampling
> > frequency Fs >= 4*f0.
>
> So, say I'm interested in frequencies between [0 and 0.1] HZ,
> my sampling frequency has to be around 0.025??

No, the maximum frequency the fft can accurately represent
is Fs/2. Therefore Fs should not be less than 2*f0. However,
f0 = Fs/2 is the last bin whereas f0 = Fs/4 is the middle bin.

Remember that Fs/f0 = T0/dt = number of samples per period.
Fs ~ 4*f0 only picks up 4 points per period.

If you had to sketch a shifted sine wave by hand in front of
an important audience, how many points per period would
you like specified?

> (say my signal is the solution of an ODE, do I have to choose t
> in a way similar to N=1000 (to say I want 1000 points),
> t=[0:1/0.025:N/0.025] ?)

t = dt*(0:N-1);
t = 0:dt:T-dt; % T = N*dt is the imposed period of ifft(fft(x))
t = linspace(0,T-dt,N);

> > 2. If abs(x[N]-x[1]) is large, I apply windowing
> > to eliminate the ficticious high frequencies
> > caused by the discontinuity
>
> How can I apply windowing? Is there any function in Matlab?
> I've an idea of what windowing does, but only read of it.

help window
doc window
lookfor window

Hope this helps.

Greg

Subject: Help with FFT

From: Luca Cerone

Date: 10 Sep, 2010 12:32:05

Message: 15 of 16

Hi Greg, I finally think I've sorted things out, but would like to be sure if my assumptions
are correct.

So, say I'm interested in knowing how frequencies in the range [0 1] behaves,
and say that I want to know them with a "resolution" of 0.01 (that is, I'm interested
in frequencies [0:0.01:1].

What I do is:
1) evaluate the frequency sampling rate
that in my case is:
fs=2; %(the maximum frequency I'm interested is one).
dt=1/fs;
Now, to be sure that I get coefficients for all the parameters
I'm interested in I set:
fr=.01 %my frequency "resolution"
and evaluate the number of samples needed:
N=fs/fr;
create the vector of times:
t=dt*[0:(N-1)];
And period of the signal is:
T=N*dt;

Now I evaluate the signal:
s=feval(myfun,t);

f=fft(s);
p=abs(f)/(N/2);
 And now, assuming N is even
freq=[0:(N/2)]/T;
p=p(1:(N/2 +1)).^2;

Does it seems correct?
Thanks again for your help,
Cheers, Luca

Subject: Help with FFT

From: Greg Heath

Date: 11 Sep, 2010 10:47:11

Message: 16 of 16

On Sep 10, 8:32 am, "Luca Cerone" <luca_cerone#_remove_th...
> @yahoo.it> wrote:
> Hi Greg, I finally think I've sorted things out, but would
> like to be sure if my assumptions > are correct.
>
> So, say I'm interested in knowing how frequencies in the
> range [0 1] behaves,

(Fs/N)*ceil((N-1)/2) = 1

> and say that I want to know them with a "resolution" of 0.01
> (that is, I'm interested in frequencies [0:0.01:1].

Fs/N = 0.01

Therefore

N = 201 % N is odd
Fs = 2.01
dt = 1/Fs
t = dt*(0:N-1);

> What I do is:
> 1) evaluate the frequency sampling rate
> that in my case is:
> fs=2; %(the maximum frequency I'm interested is one).
> dt=1/fs;
> Now, to be sure that I get coefficients for all the parameters
> I'm interested in I set:
> fr=.01 %my frequency "resolution"
> and evaluate the number of samples needed:
> N=fs/fr;
> create the vector of times:
> t=dt*[0:(N-1)];
> And period of the signal is:
> T=N*dt;
>
> Now I evaluate the signal:
> s=feval(myfun,t);
>
> f=fft(s);

Misleading notation. Usually "f" is reserved
for frequency.

Try

S = fft(s);

> p=abs(f)/(N/2);

???
Where did you get that normalization?


> And now, assuming N is even

No.

> freq=[0:(N/2)]/T;
> p=p(1:(N/2 +1)).^2;
>
> Does it seems correct?

No.
See my previous posts.

Hope this helps.

Greg

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us