fft matlab, scaling amplitude problem

1 view (last 30 days)
Bob GH
Bob GH on 30 May 2013
Hi I've faced a problem in my homework coding. please help me to solve it.
I wrote a code for a problem, which i have its results. but after taking fft function from input, the amplitude of output is twice more than expected result.
what are the reasons that might be caused this?
I really appreciate it.
  1 Comment
Bob GH
Bob GH on 30 May 2013
I wrote a part of my code related to fft :
clc; clear all;
fm=40; T=10; a=0; N=32768; %2^15 Cm=1.5; h=(T-a)/N; omega=2*pi*fm/T; hh=1/h; Iamp=2.25e3; V(1)=-65;
mhf(1)=mhfinf(V(1));
mhs(1)=mhsinf(V(1));
mNap(1)=mNapinf(V(1));
EL=leak(V(1),mhf(1),mhs(1),mNap(1));
t(1)=a;
Iinp(1)=a;
for j=1:N-1
Iinp(j+1)=Iamp*sin(omega*(t(j)^2));
%%%%%%%%% 4 step method of Runge-Kutta
K1=h*fftNapH(t(j),[V(j); mhf(j); mhs(j); mNap(j)],EL,Iinp(j+1)); deffv1=K1(1,1);mhf1=K1(2,1);mhs1=K1(3,1);mNap1=K1(4,1);
K2=h*fftNapH(t(j)+(h/2),[V(j)+(deffv1/2);mhf(j)+(mhf1/2);mhs(j)+(mhs1/2);mNap(j)+(mNap1/2)],EL,Iinp(j+1)); deffv2=K2(1,1);mhf2=K2(2,1);mhs2=K2(3,1);mNap2=K2(4,1);
K3=h*fftNapH(t(j)+(h/2),[V(j)+(deffv2/2);mhf(j)+(mhf2/2);mhs(j)+(mhs2/2);mNap(j)+(mNap2/2)],EL,Iinp(j+1)); deffv3=K3(1,1);mhf3=K3(2,1);mhs3=K3(3,1);mNap3=K3(4,1);
K4=h*fftNapH(t(j)+h,[V(j)+deffv3;mhf(j)+mhf3;mhs(j)+mhs3;mNap(j)+mNap3],EL,Iinp(j+1)); deffv4=K4(1,1);mhf4=K4(2,1);mhs4=K4(3,1);mNap4=K4(4,1);
V(j+1)=V(j)+1/6*(deffv1+(2*deffv2)+(2*deffv3)+deffv4); mhf(j+1)=mhf(j)+1/6*(mhf1+2*mhf2+2*mhf3+mhf4); mhs(j+1)=mhs(j)+1/6*(mhs1+2*mhs2+2*mhs3+mhs4); mNap(j+1)=mNap(j)+1/6*(mNap1+2*mNap2+2*mNap3+mNap4);
t(j+1)=t(j)+h;
end
V2=fft(V1,N); V3=V2(1:N/2); amp=abs(V3); ampb=amp(2:length(amp)); % remove DC part f= (0:N/2-1)*hh/N;
figure(2),plot(f,ampb) xlabel('Freq') ylabel('Amplitude') axis([0 50 0 2000])

Sign in to comment.

Answers (2)

Wayne King
Wayne King on 30 May 2013
You should always show your code:
fs = 1000;
t = 0:1/fs:1-1/fs;
L = length(x);
xdft = fft(x)/L;
plot(abs(xdft))
Exactly as I expect two peaks with amplitude 0.5
Or
xdft = 2*fft(x)/L;
xdft = xdft(1:length(x)/2+1);
plot(abs(xdft))
  4 Comments
Bob GH
Bob GH on 30 May 2013
thanks for your reply Wayne King. i wrote a relevant part of my code
Bob GH
Bob GH on 30 May 2013
I wrote a part of my code related to fft :
clc; clear all;
fm=40; T=10; a=0; N=32768; %2^15 Tfft=10; Cm=1.5; h=(T-a)/N; omega=2*pi*fm/T; hh=1/h; Iamp=2.25e3; V(1)=-65;
mhf(1)=mhfinf(V(1));
mhs(1)=mhsinf(V(1));
mNap(1)=mNapinf(V(1));
EL=leak(V(1),mhf(1),mhs(1),mNap(1));
t(1)=a;
Iinp(1)=a;
for j=1:N-1
Iinp(j+1)=Iamp*sin(omega*(t(j)^2));
%%%%%%%%% 4 step method of Runge-Kutta
K1=h*fftNapH(t(j),[V(j); mhf(j); mhs(j); mNap(j)],EL,Iinp(j+1)); deffv1=K1(1,1);mhf1=K1(2,1);mhs1=K1(3,1);mNap1=K1(4,1);
K2=h*fftNapH(t(j)+(h/2),[V(j)+(deffv1/2);mhf(j)+(mhf1/2);mhs(j)+(mhs1/2);mNap(j)+(mNap1/2)],EL,Iinp(j+1)); deffv2=K2(1,1);mhf2=K2(2,1);mhs2=K2(3,1);mNap2=K2(4,1);
K3=h*fftNapH(t(j)+(h/2),[V(j)+(deffv2/2);mhf(j)+(mhf2/2);mhs(j)+(mhs2/2);mNap(j)+(mNap2/2)],EL,Iinp(j+1)); deffv3=K3(1,1);mhf3=K3(2,1);mhs3=K3(3,1);mNap3=K3(4,1);
K4=h*fftNapH(t(j)+h,[V(j)+deffv3;mhf(j)+mhf3;mhs(j)+mhs3;mNap(j)+mNap3],EL,Iinp(j+1)); deffv4=K4(1,1);mhf4=K4(2,1);mhs4=K4(3,1);mNap4=K4(4,1);
V(j+1)=V(j)+1/6*(deffv1+(2*deffv2)+(2*deffv3)+deffv4); mhf(j+1)=mhf(j)+1/6*(mhf1+2*mhf2+2*mhf3+mhf4); mhs(j+1)=mhs(j)+1/6*(mhs1+2*mhs2+2*mhs3+mhs4); mNap(j+1)=mNap(j)+1/6*(mNap1+2*mNap2+2*mNap3+mNap4);
t(j+1)=t(j)+h;
end
V2=fft(V1,N); V3=V2(1:N/2); amp=abs(V3); ampb=amp(2:length(amp)); % remove DC part f= (0:N/2-1)*hh/N;
figure(2),plot(f,ampb) xlabel('Freq') ylabel('Amplitude') axis([0 50 0 2000])

Sign in to comment.


Azzi Abdelmalek
Azzi Abdelmalek on 30 May 2013
That means that, before calculating the fft, you've made some errors, which we can't find, because you have not posted the code.
  2 Comments
Bob GH
Bob GH on 30 May 2013
thanks for your reply. actually, It has 10++ line of code and i do not know which part do you need? is there any distinguished reason for my problem? because, my result is exactly twice of expected result.
Wayne King
Wayne King on 30 May 2013
Edited: Wayne King on 30 May 2013
Does 10++ mean 12 lines? If it is a reasonable number, please post it all. Did you look at the code I posted below?

Sign in to comment.

Categories

Find more on Fourier Analysis and Filtering 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!