fft matlab, scaling amplitude problem

1 Ansicht (letzte 30 Tage)
Bob GH
Bob GH am 30 Mai 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 Kommentar
Bob GH
Bob GH am 30 Mai 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])

Melden Sie sich an, um zu kommentieren.

Antworten (2)

Wayne King
Wayne King am 30 Mai 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 Kommentare
Bob GH
Bob GH am 30 Mai 2013
thanks for your reply Wayne King. i wrote a relevant part of my code
Bob GH
Bob GH am 30 Mai 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])

Melden Sie sich an, um zu kommentieren.


Azzi Abdelmalek
Azzi Abdelmalek am 30 Mai 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 Kommentare
Bob GH
Bob GH am 30 Mai 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 am 30 Mai 2013
Bearbeitet: Wayne King am 30 Mai 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?

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Fourier Analysis and Filtering finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by