Problem with the magnitude of DC component (zeroth order) by FFT

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

Answer 1

[Walter Roberson]

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.

[Walter Roberson]

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.

[Hwang Ryan]

Thank you for yout kind answer.

I understand that because the constant added is 5 and the duration is 1, the result is the sum of the constant 5 that I add and the product of the constant 5 and the duration 1. Right?

Answer 2

[John BG]

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

[e-mail]

{%%%%%%%% % 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)');}