fft matlab, scaling amplitude problem

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

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

1 Stimme

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. i tried it, if i apply (/L) the result is absolutely different with expected result, but if i remove (/L) the result is same as mine.
Wayne King
Wayne King am 30 Mai 2013
Again, you're not showing your code (not sure why). If you are using the fft() for amplitude estimation, you want to divide by the length of the input signal.
Bob GH
Bob GH am 30 Mai 2013
thanks for your reply Wayne King. i wrote a relevant part of my code
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

0 Stimmen

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 Get Started with Signal Processing Toolbox finden Sie in Hilfe-Center und File Exchange

Gefragt:

am 30 Mai 2013

Community Treasure Hunt

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

Start Hunting!

Translated by