Posted: March 11, 2013

If you are not familiar with Fourier Analysis, the purpose of the analysis is to represent a function as the weighted sum of a collection of trigonometric functions (sines and cosines). This, in turn, allows us to extract the spectral properties of the function.

In this case, we are going to explicitly construct the Fourier series of an input signal. This will not only give us the coefficients for each sine and cosine term in the range of frequencies we are interested in, but it will also give us the phase and magnitude of those frequencies as well.

Before diving into the Modelica code, let's be clear about what is
really going on here mathematically. Assume we have some input
signal, \(u\), that we are interested in analyzing. Further assume
that the frequencies we are interested in are the first \(n\)
harmonics of a base frequency, \(F_0\) (*i.e.,* \(F_0\), \(2 F_0\),
\(3 F_0\), ..., \(n F_0\)).

What we wish to solve for are the coefficients \(a_i\) and \(b_i\) such that:

\[u = \frac{a_0}{2} + \sum_{i=1}^{n} \left[ a_i cos(2 \pi i F_0 t) + b_i sin(2 \pi i F_0 t) \right] \]

where \(F_0\) is in Hertz. So the question is, how do we compute \(a_i\) and \(b_i\)?

Fortunately, Jean-Baptiste Joseph Fourier already worked that out for us. It turns out to be relatively simple. We integrate over the period of our base frequency as follows:

\[a_i = 2 \int_{t=0}^{\frac{1}{F_0}} u\ cos(2 \pi i F_0 t) \]

and

\[b_i = 2 \int_{t=0}^{\frac{1}{F_0}} u\ sin(2 \pi i F_0 t) \]

If we manage to work out \(a_i\) and \(b_i\), then we can compute the phase, \(\phi_i\), and magnitude, \(m_i\) using the following relationships:

\[\phi_i = tan^{-1}(b_i/a_i)\]

\[m_i = \sqrt{a_i^2+b_i^2} \]

When we want to apply all this math to the behavior of a model, all we need to do is build a simple signal processing block that performs these integrals and applies these relationships.

Let's take this bit by bit. First, what is the "public interface" of the component (the parts the user has to know about)? It should include a parameter for the fundamental frequency, a way to feed a signal in for analysis and way to extract the various analysis results as output. In other words:

```
parameter Modelica.SIunits.Frequency F0 "Base frequency for analysis";
parameter Integer n "Number of harmonics";
Modelica.Blocks.Interfaces.RealInput u "Input signal";
Modelica.Blocks.Interfaces.RealOutput a0 "Signal bias";
Modelica.Blocks.Interfaces.RealOutput a[n] "Fourier coefficients for cosine terms";
Modelica.Blocks.Interfaces.RealOutput b[n] "Fourier coefficients for sine terms";
Modelica.Blocks.Interfaces.RealOutput mag[n] "Magnitude for each frequency";
Modelica.Blocks.Interfaces.RealOutput phase[n] "Phase for each frequency";
```

Internal to the model, there are a number of things we would like to
compute. So we'll create a `protected`

section for those:

```
protected
parameter Modelica.SIunits.Time dt = 1.0/F0 "Period at base frequency";
Real s[n] = {sin(2*pi*F0*i*time) for i in 1:n} "Sine waves at various frequencies";
Real c[n] = {cos(2*pi*F0*i*time) for i in 1:n} "Cosine waves at various frequencies";
Real a0i "Integral of bias term";
Real ai[n] "Integral of cosine terms";
Real bi[n] "Integral of sine terms";
Real f "Reconstructed function";
```

where `s`

is a vector of the sine functions at each of the various
frequencies we are interested in, `c`

is a vector of cosine functions
at those same frequencies. We will use `a0i`

, `ai`

and `bi`

to hold
the values of various integrals. Finally, `f`

will compute the value
of our input signal approximated as a Fourier series (delayed by one
period of our fundamental frequency.

At the start of our simulation, we initialize all our integral terms to zero:

```
initial equation
a0i = 0;
ai = zeros(n);
bi = zeros(n);
```

We also have the following continuous equations in our models to compute the various integrals, the reconstructed function, phase and magnitude:

```
equation
der(a0i) = 2*u;
der(ai) = 2*u*c;
der(bi) = 2*u*s;
f = a0/2+a*c+b*s;
mag = {sqrt(a[i]^2+b[i]^2) for i in 1:n};
phase = {atan2(b[i], a[i]) for i in 1:n};
```

These equations demonstrate some of the interesting vector related
functions in Modelica. For example, `*`

is used as both
multiplication (scalar times vector) and as an inner product (between
vectors) in expressions like `2*u*c`

. We also use the array
comprehensions feature to compute the magnitude and phase. This
allows us to write down an expression for the \(i^{th}\) element of
a vector.

The only complicated part of the model is what we do at the end of
each period of our fundamental frequency. This is represented in the
following `when`

clause:

```
equation
// ...
when sample(0,dt) then
a0 = pre(a0i);
a = pre(ai);
b = pre(bi);
reinit(a0i, 0);
for i in 1:n loop
reinit(ai[i], 0);
reinit(bi[i], 0);
end for;
end when;
```

At the end of each period of our fundamental frequency, we extract the
current value of our integral variables (`a0i`

, `ai`

and `bi`

) and
assign their values to the coefficients in the series. Then, we use
the `reinit`

function to reinitialize those integrals so they start
from zero again.

Putting this all together, our model looks like this:

```
within Sensors.SignalProcessing;
model FourierAnalysis "Compute Fourier coefficients of an input signal"
parameter Modelica.SIunits.Frequency F0 "Base frequency for analysis";
parameter Integer n "Number of harmonics";
Modelica.Blocks.Interfaces.RealInput u "Input signal";
Modelica.Blocks.Interfaces.RealOutput a0 "Signal bias";
Modelica.Blocks.Interfaces.RealOutput a[n] "Fourier coefficients for cosine terms";
Modelica.Blocks.Interfaces.RealOutput b[n] "Fourier coefficients for sine terms";
Modelica.Blocks.Interfaces.RealOutput mag[n] "Magnitude for each frequency";
Modelica.Blocks.Interfaces.RealOutput phase[n] "Phase for each frequency";
protected
import Modelica.Constants.pi;
import Modelica.Math.atan2;
parameter Modelica.SIunits.Time dt = 1.0/F0 "Period at base frequency";
Real s[n] = {sin(2*pi*F0*i*time) for i in 1:n} "Sine waves at various frequencies";
Real c[n] = {cos(2*pi*F0*i*time) for i in 1:n} "Cosine waves at various frequencies";
Real a0i "Integral of bias term";
Real ai[n] "Integral of cosine terms";
Real bi[n] "Integral of sine terms";
Real f "Reconstructed function";
initial equation
a0i = 0;
ai = zeros(n);
bi = zeros(n);
equation
der(a0i) = 2*u;
der(ai) = 2*u*c;
der(bi) = 2*u*s;
f = a0/2+a*c+b*s;
mag = {sqrt(a[i]^2+b[i]^2) for i in 1:n};
phase = {atan2(b[i], a[i]) for i in 1:n};
when sample(0,dt) then
a0 = pre(a0i);
a = pre(ai);
b = pre(bi);
reinit(a0i, 0);
for i in 1:n loop
reinit(ai[i], 0);
reinit(bi[i], 0);
end for;
end when;
end FourierAnalysis;
```

I've included this `FourierAnalysis`

block in my open source `Sensors`

package on GitHub.

Now, I would never really undertake building a model like this without testing it. Let's first consider a case where we have only a single sine wave as an input. In this case, let's make the frequency of that sine wave three times that of our fundamental frequency. In diagram form, our test looks like this:

A simple way to determine whether the model is correct is to do a
visual comparison between the input signal and the reconstructed
function, `f`

, contained within the `FourierAnalysis`

block. This
works particularly well if each period of the input signal is
identical. In the following figure, we can see a comparison between
the input and reconstructed functions. We see clearly that they are
identical.

In the figure above, we also see the values for `b[3]`

and `mag[3]`

which should both match the magnitude of the input signal, `u`

, as
expected.

Another test we can do is to add together many different sine and
cosine waves of different magnitudes into a single waveform and then
see of our `FourierAnalysis`

block can extract them back out. That is
the approach behind our second test:

Again, since this function repeats itself with every period, we can compare the input function to the reconstructed function for a quick visual confirmation:

This is an example of how a useful analytical technique like Fourier analysis can be implemented in Modelica. Such a block can then be used to perform on-the-fly analysis while simulating a Modelica model.

As always, the source code from these examples can be found on GitHub. The code is distributed under a Creative Commons license so you can feel free to use the code in your own projects.

## Share your thoughts

comments powered by Disqus