Fourier series fit to a periodic function with multiple frequencies

Hi everyone,

I want to fit a Fourier series to a square wave functoin with two frequencies and varying amplitude. I have attached the graph for more clarification.

I can estimate the simple square wave with Fourier series. But the problem arises when I modify the amplitude of the wave with some other frequency. In the figure below, the square wave period is the same as blue function. But at some points the amplitudes are changing and the period of this part is different than the main period. Fourier series provides the average as DC part and cannot model this.

I would be grateful if you share with me your ideas on how to solve the problem.

Best,

Behrang

[Paul]

Hi Behrang,

As always, it's better to post some code and data that illustrates the problem.

Also, can you clarify the problem itself? It sounds like you wish to find the Fourier series of a signal that is expressed as the sum of two other signals. Is that correct?

Can the signal (i guess the blue curve?) be expressed as the sum of two, periodic square waves?

If so, have you verified that that signal (blue) is also periodic?

[Behrang Hoseini]

Hi Paul,

Thank you for your answer.

Well, actually this is a wave with some defect. Let's say the there is a square wave that you can easily estimate its Fourier Series counterpart. We can say that there are two waves but am not sure if can be expressed as a summation.

To generate the function, I have written the following code. it may clarify more.

I have defined a periodic window to multiply with the main wave but the problem is that the final function becomes discontinuous.

fs = 5e05;%[Hz]
f_sun =  27.4286;% [Hz]
T = 1/f_sun;
n  = 4*fix(T*fs);% number of steps   TWO SUN revolutions
Fm = 600;% [Hz]
tspan  = (0:n-1)/fs;% Solution time [s]
t    = tspan;
kl = 3e08;   %[N/m]  
ku = 5e08;  
m_is = 1.45; 
m_ir = 1.64; 
t_sun = 0.0156;

k_mean = ku*(m_is-1) + kl*(2-m_is);
k_t1 = 1/pi*(sin(2*pi*1*m_is)*cos(1*2*pi*Fm*t)+(1-cos(2*pi*1*m_is))*sin(1*2*pi*Fm*t))+...
    1/(2*pi)*(sin(2*pi*2*m_is)*cos(2*2*pi*Fm*t)+(1-cos(2*pi*2*m_is))*sin(2*2*pi*Fm*t))+...
    1/(3*pi)*(sin(2*pi*3*m_is)*cos(3*2*pi*Fm*t)+(1-cos(2*pi*3*m_is))*sin(3*2*pi*Fm*t))+...
    1/(4*pi)*(sin(2*pi*4*m_is)*cos(4*2*pi*Fm*t)+(1-cos(2*pi*4*m_is))*sin(4*2*pi*Fm*t))+...
    1/(5*pi)*(sin(2*pi*5*m_is)*cos(5*2*pi*Fm*t)+(1-cos(2*pi*5*m_is))*sin(5*2*pi*Fm*t))+...
    1/(6*pi)*(sin(2*pi*6*m_is)*cos(6*2*pi*Fm*t)+(1-cos(2*pi*6*m_is))*sin(6*2*pi*Fm*t))+...
    1/(7*pi)*(sin(2*pi*7*m_is)*cos(7*2*pi*Fm*t)+(1-cos(2*pi*7*m_is))*sin(7*2*pi*Fm*t))+...
    1/(8*pi)*(sin(2*pi*8*m_is)*cos(8*2*pi*Fm*t)+(1-cos(2*pi*8*m_is))*sin(8*2*pi*Fm*t))+...
    1/(9*pi)*(sin(2*pi*9*m_is)*cos(9*2*pi*Fm*t)+(1-cos(2*pi*9*m_is))*sin(9*2*pi*Fm*t))+...
    1/(10*pi)*(sin(2*pi*10*m_is)*cos(10*2*pi*Fm*t)+(1-cos(2*pi*10*m_is))*sin(10*2*pi*Fm*t))+...
    1/(11*pi)*(sin(2*pi*11*m_is)*cos(11*2*pi*Fm*t)+(1-cos(2*pi*11*m_is))*sin(11*2*pi*Fm*t))+...
    1/(12*pi)*(sin(2*pi*12*m_is)*cos(12*2*pi*Fm*t)+(1-cos(2*pi*12*m_is))*sin(12*2*pi*Fm*t))+...
    1/(13*pi)*(sin(2*pi*13*m_is)*cos(13*2*pi*Fm*t)+(1-cos(2*pi*13*m_is))*sin(14*2*pi*Fm*t))+...
    1/(14*pi)*(sin(2*pi*14*m_is)*cos(14*2*pi*Fm*t)+(1-cos(2*pi*14*m_is))*sin(14*2*pi*Fm*t))+...
    1/(15*pi)*(sin(2*pi*15*m_is)*cos(15*2*pi*Fm*t)+(1-cos(2*pi*15*m_is))*sin(15*2*pi*Fm*t));
k_t = (ku-kl)*k_t1 + k_mean;


i   = 1    ;% number of the planet gear, first i =  0
Tm  = 1/Fm ;
Tp  = round(3*t_sun,2);% [s] Period of the wave based on sun rotation
t1  = 0.0 + i*t_sun;
t2  = (m_is-1)*Tm + i*t_sun;% [s] mesh period of first double tooth pair  
t3  =  m_is*Tm + i*t_sun;% [s] mesh period of single tooth pair + second double tooth pair

y   = ( mod(t,Tp)>t1& round(mod(t,Tp),4)<t2)*0.95 + (round(mod(t,Tp),4)>t2 & round(mod(t,Tp),4)<t3)*0.84;
y1  = ~y;
w = y+y1;
k_t = k_t .* w;% Multiply with window W
if w <= 0
    disp('window function is zero')
    pause (5)
end 

plot(t,k_t,'r')
xlabel('t [s]');ylabel('k_m [N/m]')
set(gca,'FontSize',24) 
xlim([0 0.05])
ylim([2 6]*1e08)
grid

 

Answer 1

[Star Strider]

Adding a second function will change the baseline, for example —

t = linspace(0, 10, 1000);
s = sum((1./[1;3;5]).*sin([1;3;5]*2*pi*t*0.5)) + 0.1*cos(2*pi*t*0.1);

figure
plot(t,s)
grid

It might be appropriate to use some sort of stepwise funciton rather than a continuous function, depending on the result you want.

.

[Behrang Hoseini]

Hi Star Strider,

Thank you for your code. I have attached my own code. I got your idea. But how can I adapt it to my own code? Do you have any idea?

[Star Strider]

You apparently want to add two step plots together, with varying step amplitudes. You could probably do something like that with this example (using a random stepped baseline) —

t = linspace(0, 10, 1000);                                  % Time Vector
s = sum((1./[1;3;5]).*sin([1;3;5]*2*pi*t*0.5)) + 1;         % Synthetic Square Wave

r = max(t);                                                 % Rows Of 'ones' Matrix
c = fix(numel(t)/max(t));                                   % Columns Of 'ones' Matrix
baseline = ones(r, c) .* randi(99,r,1)/50;                  % Create 'ones' Matrix
rbaselinev = reshape(baseline.',1,[]);                      % Use 'reshape' TO Create Appropriaate Size Vector
sb = s + rbaselinev;                                        % Add Baseline Vector To Signal

figure
plot(t,s,'-b')
hold on
plot(t,sb,'-g')
plot(t,rbaselinev,'--r')
hold off
grid
legend('Original Signal','Signal With Stepped Baseline','Stepped Baseline', 'Location','best')

It would likely be relatively straightforward to adapt this approach to your problem. (I am using my code here because I understand it and so understand how it works.)

.

[Behrang Hoseini]

Thank you again.

Well, It does not give the result that I want.

Let me clarify it in another way.

I want to estimate a square wave with varying amplitude with the Fourier series.

So based on the graph below:

I. the blue function is the one that I want to estimate. There is a main period, Tm

and

II. variation of amplitude occurs at a period of Tf and it lasts for k*Tm seconds.

[Star Strider]

Using my code (for the same reasons as previously) use the fft function if you already know the stepwise baseline vector —

t = linspace(0, 10, 1000);                                  % Time Vector
s = sum((1./[1;3;5]).*sin([1;3;5]*2*pi*t*0.5)) + 1;         % Synthetic Square Wave

r = max(t);                                                 % Rows Of 'ones' Matrix
c = fix(numel(t)/max(t));                                   % Columns Of 'ones' Matrix
baseline = ones(r, c) .* randi(99,r,1)/50;                  % Create 'ones' Matrix
rbaselinev = reshape(baseline.',1,[]);                      % Use 'reshape' TO Create Appropriaate Size Vector
sb = s + rbaselinev;                                        % Add Baseline Vector To Signal

L = numel(t);
Fs = 1/(t(2)-t(1));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTbl = fft(rbaselinev,NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, mag2db(abs(FTbl(Iv))*2))
Ax = gca;
Ax.XScale = 'log';
grid
xlabel('Frequency')
ylabel('Magnitude (dB)')
title('One-Sided Fourier Transform of ‘baseline’ Signal')

figure
plot(t,s,'-b')
hold on
plot(t,sb,'-g')
plot(t,rbaselinev,'--r')
hold off
grid
legend('Original Signal','Signal With Stepped Baseline','Stepped Baseline', 'Location','best')

The Fourier transform if the ‘rbaselinev’ signal is extremely irregular because of all the abrupt transitions.

It is unlikely that any finite Fourier transform will be able to model it with any accuracy.

[Behrang Hoseini]

Thank you so much, Star Strider for your code.

Well, still it is not what I am looking for. You know I will use the wave function in an ODE for integration. Therefore, it must be smooth. As you can see in the "singal with stepped baseline", there sudden jumps where it is supposed to be smooth.

So, basically the purpose is to estimate the blue sqaure wave I shared at the beginning.

[Star Strider]

It is just not going to be ‘smooth’ in any sense with a finite Fourier transform. The best you can hope for is to window it (to compensate for the Fourier transform not being infinite) and go with that result —

t = linspace(0, 10, 1000);                                  % Time Vector
s = sum((1./[1;3;5]).*sin([1;3;5]*2*pi*t*0.5)) + 1;         % Synthetic Square Wave

r = max(t);                                                 % Rows Of 'ones' Matrix
c = fix(numel(t)/max(t));                                   % Columns Of 'ones' Matrix
baseline = ones(r, c) .* randi(99,r,1)/50;                  % Create 'ones' Matrix
rbaselinev = reshape(baseline.',1,[]);                      % Use 'reshape' TO Create Appropriaate Size Vector
sb = s + rbaselinev;                                        % Add Baseline Vector To Signal

L = numel(t);
Fs = 1/(t(2)-t(1));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTbl = fft(hann(L).*rbaselinev,NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, mag2db(abs(FTbl(Iv))*2))
Ax = gca;
Ax.XScale = 'log';
grid
xlabel('Frequency')
ylabel('Magnitude (dB)')
title('One-Sided Fourier Transform (Hanning Window) of ‘baseline’ Signal')

There are several options of the Windows. Use the one that gives the result you want.

.

[Behrang Hoseini]

Hi Star Rider,

thank you so much again for the code. Sorry I was late coming back to you.

There is a problem with this approach. You know I am going to use this function as a coefficient in an ODE. Therefore, any jumps will cause problems in the integration process. Such a jum is visible around t =8 in your code I am attaching below

So I was thinking to apply the changes on the square wave before estimating it with sine and cosine functions. But the problem is that such a amplitude varying can not be estimated with sine and cosines because of varying mean.

t = linspace(0, 10, 1000);                                  % Time Vector
s = sum((1./[1;3;5]).*sin([1;3;5]*2*pi*t*0.5)) + 1;         % Synthetic Square Wave
r = max(t);                                                 % Rows Of 'ones' Matrix
c = fix(numel(t)/max(t));                                   % Columns Of 'ones' Matrix
baseline = ones(r, c) .* randi(99,r,1)/50;                  % Create 'ones' Matrix
rbaselinev = reshape(baseline.',1,[]);                      % Use 'reshape' TO Create Appropriaate Size Vector
sb = s + rbaselinev;                                        % Add Baseline Vector To Signal

plot(t,sb, 'Linewidth',1.0)
grid 

[Star Strider]

‘But the problem is that such a amplitude varying can not be estimated with sine and cosines because of varying mean.’

Absoslutely. If you are going to estimate parameters, especially if using them with differential equations, it will be necessary to use them in a stepwise approach, starting and stopping the estimation procedure at the beginning and end of each step.

[Behrang Hoseini]

For the integration I am using Generalize alpha method which is kind of a generalization of Newmark method. It uses varying time steps. So when I use stepwise functions it causes instability.

Is that possible to estimate such a function with Wavelets?

[Star Strider]

I am lost now with respect to what you are doing. This began as an entirely different problem. I am not certain what the best approach would be.

[Behrang Hoseini]

Sorry, I should have explained it in a better way.

What I want to to is to solve the following ODE

the K parameter is varting with time like a stepwise function as below

Generally it is estimated with a Fourier Series.

but in some cases, the amplitude varies with time and won't be fully periodic with a single frequency. This is the problem I am facing. I tried using a window and multiply it with the wave but due to the jumps, it causes instability in the integration.

I guess it is clear now.

[Star Strider]

I have no further suggestions. What I posted is the best I can do.

[Behrang Hoseini]

OK. Thank you so much. Actually, I learned so many things from the code you shared.

[Star Strider]

My pleasure!

Answer 2

[Benjamin Thompson]

Try looking over the documentation of the fft in MATLAB, along with the examples included there. If you have more specific questions and can post some of your work that you have a question about then post that to this board or add it to this questions.

[Behrang Hoseini]

Hi Benjamin,

Thank you for your answer. Well, I have checked it before. I couldn't get my answer. Now I have attached my own code. Can you please see it.

Best,

Behrang

[Benjamin Thompson]

If you plot the FFT result, you will see there are many harmonics needed to approximate both the square wave part and the odd change amplitude or bias. So two frequencies added together will not approximate your signal very well.

f = fs*(0:(n/2 - 1))/n;

fft_kt = abs(fft(k_t)/n);

figure, plot(f, fft_kt(1:n/2)), grid on, zoom on

[Behrang Hoseini]

Thank you Benjamin,

Well you are right. My aim is to estimate the square wave by any method. I have used 15 harmonics to estimate the wave. However, it does not work for the case of varying amplitude.