How do I use a for-loop to do fft and plot the power?

Hello,

I have a data set (time and variable Be) and I want to do a fft on one part of the data at a time (1000 years). To do this, I would like to use a for-loop that plots the power of the period for each of the different ffts. How do I do this? I'm new to Matlab but I have tried the code below:

hold on 

for ts=(0:1000:5000)';             % I have data between 0 and 5000 years ago. 
    A=fft(Be);                     % Do FFT on the Be data.
    A(1)=[];
    n=length(A);                     
    power=abs(A(1:floor(n/2))).^2;
    nyquist=1/2;
    freq=(1:n/2)/(n/2)*nyquist;
    period=1./freq;
    plot(period,power);
end

When I do this, I get only one plot. What am I doing wrong?

Thanks!

Answer 1

[Thorsten]

May be that's because they are the same? I get different curves using random input.

BTW: note that n is always 1000 in your case, so you can compute "period" outside the loop; further, as you have an offset of 350 now, a could start from 0; just be(0) is wrong syntax in Matlab because matrix indices start with 1.

figure; 

be_i_dt = rand(1, 20000);

hAxes = gca;
hold( hAxes, 'on' )

offset = 305;
n = 1000;

n=length(A);                     

nyquist=1/2;
freq=(1:n/2)/(n/2)*nyquist;
period=1./freq;

for a = 1:9
    A = fft(be_i_dt(offset+a*n:offset+(a+1)*n));                 
    A(1)=[];
    power=abs(A(1:floor(n/2))).^2;
    plot(hAxes,period,power)
end

[Emma]

Thanks! This works!!

Before, I did a fft on a matrix (time and be-data) but in your code I noticed that you put a vector in the fft. So I tried that (removed the time from my matrix) and it works almost fine now.

I get an an error message, index exceeds matrix dimensions and my plot only shows 7 of the plots (still better than the 9 of the same plots). Why is this? Do I have to few data points maybe?

[Thorsten]

Yes, you probably have to few data points. The last value

 offset+(a+1)*n

for the highest a (9 in your case) in your for loop must always be lower or equal to

numel(be_i_dt)

[Emma]

Thank you!

Answer 2

[David Young]

The variable Be is not changed between iterations, so A will also be the same each time, so each plot will be the same as the previous one. If Be is a vector with all 5000 data points, then you could replace the first line in the loop with

A = fft(Be(ts+1:ts+999));

to select 1000 data points starting from ts.

[Emma]

Thanks, that explains a lot! But I tried this and still only get one plot (it looks a bit different though). Is it possible to make one plot and have five different power-plots in it? So I can see the difference between the five periods?

[Adam]

I usually prefer to use instructions such as

figure; hAxes = gca;

hold( hAxes, 'on' )

plot( hAxes, period, power )

when dealing with axes as it can be easy to have things go wrong relying on implicit behaviour of the axes you expect being the one currently in focus. That may not be the problem here as your 'plot' instructions should apply to the same axes as your 'hold' instruction.

When you say you have only one plot do you mean you can only see 1 plot or that there is literally just 1 plot in the plot browser? (i.e. has it plotted the same thing 5 times still or has it really only plotted one?).

Incidentally your loop has 6 iterations rather than 5.

[Emma]

Thanks for the tips!

There is only one plot. I did plot the ffts for each period separately (without the loop) and they look different, so the 5 plots should not look the same.

I also tried to add:

thisfig = figure();

to try to get 5 different plots, but then I get the same one as before plus an empty figure.

How come it has 6 iterations?

[Thorsten]

Could you please show your present code, and provide the Be matrix so that we can try it and see what went wrong?

Answer 3

[Sudharsana Iyengar]

HI,

Try this and let me know.

     for a= 0:4
    A=fft(Be(a*1000:(a+1)*1000));                 
    A(1)=[];
    n=length(A);                     
    power=abs(A(1:floor(n/2))).^2;
    nyquist=1/2;
    freq=(1:n/2)/(n/2)*nyquist;
    period=1./freq;
    figure()
    plot(period,power)
    end

[Emma]

Hi!

Thanks! I tried it but got an error message (Subscript indices must either be real positive integers or log) so I changed a to

a=1:5

and then it plots 5 different figures (no error message=good) but it's the same plot in each figure. So the plot doesn't change.

I have changed some things (the period I'm lookin at is now from year 305 to 9305 before present, be is now called be_i_dt) and the code looks like this now:

for a= 1:9
    A=fft(be_i_dt(305+a*1000:305+(a+1)*1000));                 
    A(1)=[];
    n=length(A);                     
    power=abs(A(1:floor(n/2))).^2;
    nyquist=1/2;
    freq=(1:n/2)/(n/2)*nyquist;
    period=1./freq;
    figure()
    plot(period,power)
end

And I still get 9 figures that look the same.