**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Fourier series fit to a periodic function with multiple frequencies

32 views (last 30 days)

Show older comments

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

##### 2 Comments

Paul
on 12 Aug 2022

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
on 22 Aug 2022

Edited: Behrang Hoseini
on 22 Aug 2022

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

### Answers (2)

Star Strider
on 12 Aug 2022

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.

.

##### 14 Comments

Behrang Hoseini
on 22 Aug 2022

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
on 22 Aug 2022

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
on 22 Aug 2022

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
on 22 Aug 2022

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
on 23 Aug 2022

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
on 23 Aug 2022

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')

.

Behrang Hoseini
on 2 Sep 2022

Edited: Behrang Hoseini
on 2 Sep 2022

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
on 2 Sep 2022

‘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
on 2 Sep 2022

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
on 2 Sep 2022

Behrang Hoseini
on 2 Sep 2022

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.

Behrang Hoseini
on 3 Sep 2022

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

Benjamin Thompson
on 12 Aug 2022

##### 3 Comments

Behrang Hoseini
on 22 Aug 2022

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
on 22 Aug 2022

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
on 23 Aug 2022

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.

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)