On Aug 20, 1:40 pm, "shinchan " <shinchan75...@gmail.com> wrote:
> I have 300 data points collected at 1Hz. I want to do FFT on this time series, then square it to get to power spectrum. Should I pad it with zero to make it 512? fft function lets me do it either way, but the fft results are different. So I don't quite know what's going on.
I have recently done a bunch of FFTs in Matlab, and I agree the
documentation is really confusing. I went ahead and followed Numerical
Recipes and wrote my own code in 1D then figured out what Matlab was
doing, which is, well, confusing. Here's some code to get you started,
but I highly recommend you read up on it:
% My Fourier Transform
close all
clear all
% h(x) = ....
% N is number of gridpoints, resolution of h(x)
% make up data:
dx = 0.1;
x = 0:dx:50.1;
% N = 50; % has to be even
% n = N/2:N/2;
% f = n./(N*dx); % from 5 to 5 in this case
N = 50;
n = 0:N1;
f = n./(N*dx); % = n/Lx
% f = 1.6:
h = cos(2*pi*x*3) + sin(2*pi*x*1.6)+2;
% cos and constant go to real, sin goes to imag H
figure
plot(x,h,'k.')
% plot(x(1:10),h(1:10))
% hold on
% plot(x(1:10),h(1:10),'k.')
xlabel('x')
ylabel('h(x)')
H = x;
H(:) = 0;
sum = 0;
count = 0;
%for c = N/2:1:N/2
for c = 0:N1
count = count +1;
for k = 0:N1
sum = sum + h(k+1)*exp(2*pi*i*k*c/N);
end
%sum = sum*dx;
H(count) = sum;
sum =0;
end
figure;
subplot(2,1,1)
%plot(f,real(H))
plot(n,real(H))
ylabel('re(H)')
hold on
%plot(f(34),real(H(34)),'*')
subplot(2,1,2)
%plot(f,imag(H),'k')
plot(n,imag(H),'k')
xlabel('n')
ylabel('im(H)')
hold on
%plot(f(34),imag(H(34)),'k*')
% Inverse finite transform:
count = 0;
%for c = N/2:1:N/2
for c = 0:N1
count = count +1;
for k = 0:N1
sum = sum + H(k+1)*exp(2*pi*i*k*c/N);
end
%sum = sum/(N*dx);
sum = sum/N;
h_back(count) = sum;
sum = 0;
end
figure(1)
hold on
plot(x,h_back,'b.')
plot(x,imag(h_back),'g')
% reproduces the data perfectly
% now check Matlab:
H_m= fft(h);
figure;
subplot(2,1,1)
plot(n,real(H_m),'r')
ylabel('re(H_m)')
hold on
subplot(2,1,2)
plot(n,imag(H_m),'r')
xlabel('f')
ylabel('im(H_m)')
hold on
% now back w/ Matlab:
h_m = ifft(H_m);
figure
plot(x,h_m,'k.')
% plot(x(1:10),h(1:10))
% hold on
% plot(x(1:10),h(1:10),'k.')
xlabel('x')
ylabel('h_m(x)')
hold on;
plot(x,h,'ro')
figure(1)
hold on
h_back_m = ifft(H);
plot(x,h_back_m, 'gs')
