Problem with the magnitude of DC component (zeroth order) by FFT
6 views (last 30 days)
Show older comments
I coded below.
I expect to get 5 as the magnitude of zeroth order. But the magnitude of DC component appears as 10 even though the magnitude of 1st and 5th order as I expected. Please let me know what the problem is and how can I solve it.
t1=[0:0.02:0.98];
y1 = 5 + 3*sin(2*pi*t1) + sin(5*2*pi*t1);
fft1 = fft(y1);
figure(1)
plot(t1,y1)
grid on
figure(2)
plot(abs(fft1)*2/length(t1));
xlabel('Harmonic order')
ylabel('Amplitude')
grid on
0 Comments
Accepted Answer
Walter Roberson
on 22 Jan 2011
The first output element reflects the total power of the system. If you were to examine the first output element when your DC offset is 0, and then to examine the first output element when your DC offset is P, you would found that the first output element has increased by P*N where N is the length of the vector.
2 Comments
Walter Roberson
on 22 Jan 2011
Remember the first component is the integral of the signal. If the signal is symmetric around 0, and thus has no DC component, then the integral is 0 and the first output would be 0. If though you add a constant to the signal, then the integral would be increased by the constant times the duration of the signal, which is the constant times the number of points for the discrete form.
More Answers (1)
John BG
on 20 Nov 2014
Edited: John BG
on 20 Nov 2014
playing a bit with your script:
{t1=[0:0.005:12]; % total amount of time samples: 10*1/.001=1001 y1 = 5 + 3*sin(2*pi*t1) + sin(2*pi*5*t1); % f1=1 f2=5 Y1fft = fft(y1); K=40;k=-K:K;w1=pi*k/K; % total amount of frequency samples: kN=2*50+1(DC) figure(1);plot(y1);grid on figure(2);Y1dft=dtft(y1,t1,k); % modY1dft=abs(Y1dft) figure(3); plot(abs(Y1fft));grid on title('|Y1dft|') % plot(abs(fft1)*2/length(t1)); xlabel('|Y1fft|');ylabel('Amplitude');grid on}
The initial t =0:.001:10 k=-50:50 DFT
DC f1 A1 f2 A2
tN=10*1000 495.1 .9549 107.9 4.934 20.77
FFT
tN DC f1 A1 f2 A2 tN/kN=1001/101 =9.91
10*100 5001x 11 15000 51 4999
1.11 5.15
for the last set I tried t=[0:.005:12] k=[-40:40]
FFT DC A1 A2
14000 3601 1198
tN/kN=12*1/.005=2400 5.83 3601/2400=1.5004 1198/2400=.499
DC a bit of course but A1 and A2 seem ok, don't they?
regards
John
{%%%%%%%% % function [X]=dtft(x,n,w) % from Digital Signal Processing, 3rd ed, by J.Proakis V.Ingle % X: DTFT values at vector w frequencies % x: finite duration sequence over n % n: sample position vector % w: frequency location vector may come as [0,pi] or [-pi,pi] % w is actually f/fo, it may be low pass or band pass signal. % let be band limited signal x(t), BW=2*pi*f0. Sampling at fs=2f0 % DTFT f goes from -fs/2 to fs/2. w goes from -ws/2 to ws/2. % % if w proportional to [0,pi] then k=w*(k_max/pi) % where k_max=length(k)-1=length(w)-1 % if w proportional to [-pi,pi] then k=w*(k_max/(2*pi)) % where k_max=(length(w)-1)/2=(length(k)-1)/2 % % input vector w should be like [0 fstep 2*fstep .. wN-fstep wN] % that means or % [-wN -wN+f_step -wN+2+2*fstep .. -fstep 0 fstep .. wN-fstep wN] % function [X,handle_graphs_X]=dtft(x,n,w)
L_w=length(w); k=zeros(1,length(w));
if ~k(1), k=((L_w-1)/pi)*w;;K2_w=length(w)-1; % there is no k=0, means w [0,w0] else k=[-((L_w-1)/(2*pi)):1:((L_w-1)/(2*pi))].*w; K2_w=(length(w)-1)/2; % w [-w0,w0] end
X=x*(exp(-j*pi/K2_w)).^(n'*k);
% X=W*x % W=[exp(-j*pi/M).*(k'*n)] % % X'=x'*W'=x'*[exp(-j*pi/M).^(n'*k)] %
%%%%% the following is verification only magX=abs(X);argX=angle(X);ReX=real(X);ImX=imag(X);
% handle_graphs_X=figure(newplot); handle_graphs_X=newplot; % hold all; subplot(2,2,1);plot(w/(2*pi),magX./length(w));grid on; hold all; xlabel('w/(2*pi)');ylabel('mod(X)'); subplot(2,2,2);plot(w/(2*pi),argX);grid on; hold all; xlabel('w/(2*pi)');ylabel('ang(X)'); subplot(2,2,3);plot(w/(2*pi),ReX);grid on; hold all; xlabel('w/(2*pi)');ylabel('Re(X)'); subplot(2,2,4);plot(w/(2*pi),ImX);grid on; hold all; xlabel('w/(2*pi)');ylabel('Im(X)');}
0 Comments
Communities
More Answers in the Power Electronics Control
See Also
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!