How do I use a for-loop to do fft and plot the power?
2 views (last 30 days)
Show older comments
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!
0 Comments
Accepted Answer
Thorsten
on 17 Dec 2014
Edited: Thorsten
on 17 Dec 2014
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
3 Comments
Thorsten
on 17 Dec 2014
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)
More Answers (2)
David Young
on 16 Dec 2014
Edited: David Young
on 16 Dec 2014
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.
4 Comments
Sudharsana Iyengar
on 17 Dec 2014
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
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!