# Fourier series fit to a periodic function with multiple frequencies

36 views (last 30 days)
Behrang Hoseini on 12 Aug 2022
Commented: Star Strider on 3 Sep 2022
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 CommentsShowHide 1 older comment
Behrang Hoseini on 22 Aug 2022
Edited: Behrang Hoseini on 22 Aug 2022
Hi Paul,
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

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.
.
Star Strider on 3 Sep 2022
My pleasure!

Benjamin Thompson on 12 Aug 2022
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 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.

### Categories

Find more on Descriptive Statistics in Help Center and File Exchange

R2021b

### Community Treasure Hunt

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

Start Hunting!

Translated by