Reconstruct a signal by determining the N Fourier Coefficients

[MoonDiver]

Homework Statement

I want to determine by using matlab the Fourier Series Complex Coefficients of the signal x(t), based on the below formula and to plot the amplitude spectre. Also, I want to rebuild the initial signal using the determined Fourier Coefficients, and I want to plot on the same graph the initial signal and the rebuilt one so that I can compare them. Unfortunately, haven't got the expected results.

Relevant Equations

x(t)=1/P*[Σ][/k][/∞]=-∞[X][/k]*[e][/jkw0t] with w0=2π/P;

[X][/k]=[∫][/P] x(t)*[e][/-jkw0t] dt for k coefficients;

[x][/^](t)=1/P*(X0+2*[Σ][/N][/k=1] [X][/k]*[e][/jkw0t] ); [X][/-k]=[X][/k][/*]; for reconstruction of the signal.

%My code:
%Type of signal: square

T = 40; %Period of the signal [s]

F=1/T;   % fr

D = 23; % length of signal(duration)
dt=(D/T)*100;
N = 50; %Number of coefficientsw0 = 2*pi/T; %signal pulset1= 0:0.002:T; % original signal sampling

x1 = square((2*pi*F)*(t1),dt);%initial square signal

t2= 0:0.002:D; %modified signal sampling

x2 = zeros(1,length(t2)); %initializing the modified signal with null values.
dif=T-D;

x2(t1<=D)=x1(t1<=D);% modify the null values with values from the original signal.
x2(1,dif:D)=x1(1,dif:D); %modify for values of t1>=T-D.subplot(2,1,1)
plot(t2,x2),title('x(t)+ reconstructed signal)');
hold on

for k = -N:N %k represents the variable after which the sum is achieved

   
x3 = x1; %x3 represents the signal obtained after the Fourier Series formula;

   
x3 = x3 .* exp(-1i*k*w0*t1);
   
X(k+N+1) = 0; %initialise with null value

end

   
for i = 1:length(t1)-1   
X(k+N+1) = X(k+N+1) + (t1(i+1)-t1(i)) * (x3(i)+x3(i+1))/2; %reconstruction using the coefficients   

endfor i = 1:length(t1)

   
x_rec(i) = 0; %initialise with null value% x_rec is the reconstructed signal using N Coefficients
end

for k=-N:N

       
x_rec(i) = x_rec(i) + (1/T) * X(k+51) * exp(1i*k*w0*t1(i));  %reconstruction using the coefficients ( the integral being calculated as a sum)
end

plot( t1, x_rec, '--')
subplot(2,1,2)
w=-50*w0:w0:50*w0; %w is the vector which allows displaying the spectre of the function

stem(w/(2*pi),abs(X));

My result:

 

[LeafNinja]

I got lost as to what you were doing around line 26 but if what you're doing is obtaining Fourier series coefficients and then reconstructing a wave then the following should help.

Continuous function Fourier series coefficients is given by
$a_{n} = \frac{1}{T} \int_0^T X(t)e^{j \frac{-2 \pi nt}{T}}dt$

In discretized form this is aprroximately equal to
$a_{n} = \frac{1}{N} \sum_{k = 0}^{N-1} X(k)e^{j \frac{-2 \pi nk}{N}}$

where k is an integer that represents time t and N is an integer that represents Period T

Similarly, when reconstructing the wave using sin waves, we have

$X(t) = \sum_{n = 0}^{N-1} a_{n} e^{j \frac{2 \pi nt}{T}}$

And the discrete form is

$X(k) = \sum_{n = 0}^{N-1} a_{n} e^{j \frac{2 \pi nk}{N}}$

[Mentor Note -- OP has made corrections to the above equations in a later post]

I would like to point out that the above formulas contain slight errors and will make the corrections shown below.

The Continuous Fourier series coefficients are given by
$a_{n} = \frac{1}{T} \int_0^T X(t)e^{j \frac{-2 \pi nt}{T}}dt$ for $-N<= n<=N$

And to reconstruct the continuous signal
$X(t) = \sum_{n = -N}^{N} a_{n} e^{j \frac{2 \pi nt}{T}}$

The following plots demonstrate the corresponding code shown below. The plots are original wave, amplitudes of different frequencies, and reconstructed wave respectively.

[CODE lang="matlab" title="Code for Fourier Series"]T=40; %Period
t = (0:.1:T); %Time vector

dt =23/40 * 100; %Duty Cycle
X1 = square( t*2*pi/T, dt); % Initial Wave

subplot(3, 1, 1);
plot(t, X1);

%Obtain Fourier cofficients

N = length(t); % Number of cofficients

k = (0:N-1); % Discretized sample time

A = zeros(1,N);
for n=0:N-1
A(n+1) = 1/N * sum(X1 .* exp(-1j*2*pi*n *k/N )); %Obtain nth amplitude coefficient
end

f = (0:N-1)* 1/T; %convert coefficient number to frequency

%Plot amplitude vs frequency
subplot(3,1,2)
absA = abs(A);
bar(f(1:int32(N/2)), absA(1:int32(N/2))); % high frequency amplitudes are usually not accurate for Fourier series because sampling frequency approaches size of Fourier series frequency
xlabel("Frequency");
ylabel("Amplitude");

%Reconstruct Original wave using Fourier series

X2 = zeros (1, N);
n = (0:N-1); %n is now a vector representing all n between 0 and N
for i=1:N
k = i - 1; %k is now discrectized time scalar
X2(i) = sum(real(A .* exp(1j*2*pi*n *k/N))); %Compute amplitude at each point in vector t,make sure it is real
end

%Plot reconstructed wave

subplot(3, 1, 3);
plot(t, X2);
[/CODE]

 

[MoonDiver]

I got lost as to what you were doing around line 26 but if what you're doing is obtaining Fourier series coefficients and then reconstructing a wave then the following should help.

Continuous function Fourier series coefficients is given by
$a_{n} = \frac{1}{T} \int_0^T X(t)e^{j \frac{-2 \pi nt}{T}}dt$

In discretized form this is aprroximately equal to
$a_{n} = \frac{1}{N} \sum_{k = 0}^{N-1} X(k)e^{j \frac{-2 \pi nk}{N}}$

where k is an integer that represents time t and N is an integer that represents Period T

Similarly, when reconstructing the wave using sin waves, we have

$X(t) = \sum_{n = 0}^{N-1} a_{n} e^{j \frac{2 \pi nt}{T}}$

And the discrete form is

$X(k) = \sum_{n = 0}^{N-1} a_{n} e^{j \frac{2 \pi nk}{N}}$

The following plots demonstrate the corresponding code shown below. The plots are original wave, amplitudes of different frequencies, and reconstructed wave respectively.
View attachment 254756
[CODE lang="matlab" title="Code for Fourier Series"]
T=40; %Period
t = (0:.1:T); %Time vector

dt =23/40 * 100; %Duty Cycle
X1 = square( t*2*pi/T, dt); % Initial Wave

subplot(3, 1, 1);
plot(t, X1);

%Obtain Fourier cofficients

N = length(t); % Number of cofficients

k = (0:N-1); % Discretized sample time

A = zeros(1,N);
for n=0:N-1
A(n+1) = 1/N * sum(X1 .* exp(-1j*2*pi*n *k/N )); %Obtain nth amplitude coefficient
end

f = (0:N-1)* 1/T; %convert coefficient number to frequency

%Plot amplitude vs frequency
subplot(3,1,2)
absA = abs(A);
bar(f(1:int32(N/2)), absA(1:int32(N/2))); % high frequency amplitudes are usually not accurate for Fourier series because sampling frequency approaches size of Fourier series frequency
xlabel("Frequency");
ylabel("Amplitude");

%Reconstruct Original wave using Fourier series

X2 = zeros (1, N);
n = (0:N-1); %n is now a vector representing all n between 0 and N
for i=1:N
k = i - 1; %k is now discrectized time scalar
X2(i) = sum(real(A .* exp(1j*2*pi*n *k/N))); %Compute amplitude at each point in vector t,make sure it is real
end

%Plot reconstructed wave

subplot(3, 1, 3);
plot(t, X2);
[/CODE]

Hi,LeafNinja ! Thank you for the response! I'm trying to reconstruct the signal but using only a limited number of coefficients. In my case N = 50.
Should obtain something like the example in the photo.

 

[LeafNinja]

Looking at the plot on the right, I cannot tell what N is but it appears that -N<=n<=N but only -5<=n<=5 is shown.

You can find the nth coefficient using the discrete Fourier transform using the equation below and then taking the absolute value. $a_{n} = \sum_{k = -N}^{N} X(k) e^{-\frac{ 2j\pi nk}{2N}}$

where -N<=k<=N corresponds to -T/2<t<T/2 and the period T appears to be 2 in this case, looking at the plot on the left.

*This is a slight variation to the discrete Fourier series formula shown above

Then, by taking the inverse discrete Fourier transform, you can reconstruct the signal. Remember that because you are using positive and negative n, the coefficients are conjugates so you only need either positive or negative.

 

[MoonDiver]

Looking at the plot on the right, I cannot tell what N is but it appears that -N<=n<=N but only -5<=n<=5 is shown.

You can find the nth coefficient using the discrete Fourier transform using the equation below and then taking the absolute value.$a_{n} = \sum_{k = -N}^{N} X(k) e^{-\frac{ 2j\pi nk}{2N}}$

where -N<=k<=N corresponds to -T/2<t<T/2 and the period T appears to be 2 in this case, looking at the plot on the left.

*This is a slight variation to the discrete Fourier series formula shown above

Then, by taking the inverse discrete Fourier transform, you can reconstruct the signal. Remember that because you are using positive and negative n, the coefficients are conjugates so you only need either positive or negative.

Thank you! The problem is that I don t understand how should I compute the Amplitude of N coefficients in my case 50, if the length of t would be equal to -T/2:0.1:T/2 and the length of k in this case is -50:50. This gives me an error because the signal x would have the length of t and I would need to multiply it with the formula described by you, but it has the length of k instead of t.

 

[MoonDiver]

Thank you again, LeafNinja! Got it now!

 

[LeafNinja]

I got lost as to what you were doing around line 26 but if what you're doing is obtaining Fourier series coefficients and then reconstructing a wave then the following should help.

Continuous function Fourier series coefficients is given by
$a_{n} = \frac{1}{T} \int_0^T X(t)e^{j \frac{-2 \pi nt}{T}}dt$

Similarly, when reconstructing the wave using sin waves, we have
$X(t) = \sum_{n = 0}^{N-1} a_{n} e^{j \frac{2 \pi nt}{T}}$
...

I would like to point out that the above formulas contain slight errors and will make the corrections shown below.

The Continuous Fourier series coefficients are given by
$a_{n} = \frac{1}{T} \int_0^T X(t)e^{j \frac{-2 \pi nt}{T}}dt$ for $-N<= n<=N$

And to reconstruct the continuous signal
$X(t) = \sum_{n = -N}^{N} a_{n} e^{j \frac{2 \pi nt}{T}}$