fft matlab, scaling amplitude problem

1 次查看（过去 30 天）
Bob GH
Bob GH 2013-5-30
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 个评论
Bob GH
Bob GH 2013-5-30
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])

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Wayne King
Wayne King 2013-5-30
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 个评论
显示 2更早的评论
Bob GH
Bob GH 2013-5-30
thanks for your reply Wayne King. i wrote a relevant part of my code
Bob GH
Bob GH 2013-5-30
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])

请先登录，再进行评论。

Azzi Abdelmalek
Azzi Abdelmalek 2013-5-30
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 个评论
Bob GH
Bob GH 2013-5-30
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 2013-5-30
编辑：Wayne King 2013-5-30
Does 10++ mean 12 lines? If it is a reasonable number, please post it all. Did you look at the code I posted below?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by