Reproducing signal by it's first 40 harmonics

I have a following periodic signal:

s(t) = {E*e^(t-ti), t <= 0; 0, t > 0}; E = 10, ti = 1; T = 10

I built a plot with 3 periods of the signal like so:

E = 10; ti = 1; T = 10;
fs = 100;
t1 = -2*T:1/fs:-T; t2 = -T:1/fs:0; t3 = 0:1/fs:T;
s1 = E*exp(t2/ti);
t = [t1 t2 t3]; s = [s1 s1 s1];

figure();
plot(t,s);

I need to find first 40 harmonics of this signal, and restore the signal by them.

How can I do that?

Answer 1

[Mathieu NOE]

hello

try this

I computed first your signal frequency then the ahrmonics complex amplitudes . From there you can easily generate the 40 harmonics truncated signal

E = 10; ti = 1; T = 10;
fs = 100;
t1 = -2*T:1/fs:-T; t2 = -T:1/fs:0; t3 = 0:1/fs:T;
s1 = E*exp(t2/ti);
t = [t1 t2 t3]; s = [s1 s1 s1];

% signal crossing point : compute signal period
threshold = max(s)/2; % 
t0_pos1 = find_zc(t,s,threshold);
period = mean(diff(t0_pos1)); % delta time of crossing points
freq = 1./period; % signal frequency

% compute dc / real (cos) and imaginary (sin) harmonics
samples = numel(t);
dc_value = mean(s);
harmo = 40;

for k = 1:harmo
    rea(k) = 2/samples*trapz(s.*cos(2*pi*k*freq*t));
    ima(k) = 2/samples*trapz(s.*sin(2*pi*k*freq*t));
end
    
% re-create signal from dc + harmonics
out = dc_value; % first load the dc value
for k = 1:harmo % then add the  harmonics
    out = out + rea(k).*cos(2*pi*k*freq*t) + ima(k).*sin(2*pi*k*freq*t);
end

figure(1)
plot(t,s,'b',t,out,'m',t0_pos1,threshold*ones(size(t0_pos1)),'*r','linewidth',2,'markersize',12);grid on
title('Signal plot');
legend('input signal','signal reconstructed','crossing points');
xlabel('time (s)')
ylabel('Amplitude')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Nikita]

edit: period and freq are NaN

Matrix dimensions must agree.

rea(k) = 2 / samples * trapz(s.*cos(2*pi*k*freq*t));

[Mathieu NOE]

funny

I simply used your code and expanded it ...

you can also prove that it works as anybody can execute my code on this web page (click the green arrow RUN)

this is how you get the plot above with the magenta curve

did you copy my entire code , with the function find_zc as it appears at the end ?

E = 10; ti = 1; T = 10;
fs = 100;
t1 = -2*T:1/fs:-T; t2 = -T:1/fs:0; t3 = 0:1/fs:T;
s1 = E*exp(t2/ti);
t = [t1 t2 t3]; s = [s1 s1 s1];
% signal crossing point : compute signal period
threshold = max(s)/2; % 
t0_pos1 = find_zc(t,s,threshold);
period = mean(diff(t0_pos1)); % delta time of crossing points
freq = 1./period; % signal frequency
% compute dc / real (cos) and imaginary (sin) harmonics
samples = numel(t);
dc_value = mean(s);
harmo = 40;
for k = 1:harmo
    rea(k) = 2/samples*trapz(s.*cos(2*pi*k*freq*t));
    ima(k) = 2/samples*trapz(s.*sin(2*pi*k*freq*t));
end
    
% re-create signal from dc + harmonics
out = dc_value; % first load the dc value
for k = 1:harmo % then add the  harmonics
    out = out + rea(k).*cos(2*pi*k*freq*t) + ima(k).*sin(2*pi*k*freq*t);
end
figure(1)
plot(t,s,'b',t,out,'m',t0_pos1,threshold*ones(size(t0_pos1)),'*r','linewidth',2,'markersize',12);grid on
title('Signal plot');
legend('input signal','signal reconstructed','crossing points');
xlabel('time (s)')
ylabel('Amplitude')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end

[Nikita]

of course I copied entire code, including find_zc function

[Mathieu NOE]

ok , so how can it be that the code works on this forum page and not on your pc ?

[Nikita]

Maybe because I'm using R2020b?

[Walter Roberson]

@Nikita

I ran the code from https://www.mathworks.com/matlabcentral/answers/2040456-reproducing-signal-by-it-s-first-40-harmonics#comment_2942146 on my Mac, R2020b, and got the same plot that @Mathieu NOE shows.

Please show the complete error message you get. Also after you get the error message please show the output of

whos

Answer 2

[Nikita]

Here's an error log:

Matrix dimensions must agree.

Error in Untitled (line 35)
    rea(k) = 2 / samples * trapz(s.*cos(2*pi*k*freq*t));

And here's whos output:

      Name           Size              Bytes  Class     Attributes

  E              1x1                   8  double              
  N              1x1                   8  double              
  N1             1x1                   8  double              
  S              1x512              8192  double    complex   
  S_mag          1x512              4096  double              
  S_pha          1x512              4096  double              
  T              1x1                   8  double              
  dc_value       1x1                   8  double              
  f              1x512              4096  double              
  freq           1x1                   8  double              
  fs             1x1                   8  double              
  harmo          1x1                   8  double              
  k              1x1                   8  double              
  period         1x1                   8  double              
  s              1x1001             8008  double              
  s1             1x1001             8008  double              
  samples        1x1                   8  double              
  t              1x3003            24024  double              
  t0_pos1        1x1                   8  double              
  t1             1x1001             8008  double              
  t2             1x1001             8008  double              
  t3             1x1001             8008  double              
  threshold      1x1                   8  double              
  ti             1x1                   8  double        

[Walter Roberson]

The code has

s = [s1 s1 s1];

Because s1 is repeated 3 times, the final result has to have a number of columns that is divisible by 3.

 s              1x1001             8008  double  

But 1001 is not divisible by 3. Therefor the code you are executing must not be the same as the one posted by @Mathieu NOE

Put a breakpoint in the code at the line

 t = [t1 t2 t3]; s = [s1 s1 s1];

and run to there, and check size(s1) and then dbstep to execute the line and then check size(s) . If size(s) does not end up exactly 3 times the size of s1 then you have a serious problem with your MATLAB. If s1 ends up length 1001 and s ends up length 3003 (same as t) then use dbstep repeatedly stepping one line at a time, checking size(t) and size(s) each line, looking for the place that s changes size down to 1001 even though there are no further assignments to s .

[Nikita]

yeah, i found the problem, fixed it, and it works, thank you