-
Coefficients of Fourier series
Ching-Tang
Tseng
Hamilton, New
Zealand
2 May 2018
Artificial
music output in my computer was synthesized audio. Off the shelf
software package information shows me that sound was mere programmed
by combination of sine and cosine functions. It means that all sound
waveform of musical instruments could be synthesized by sine or
cosine functions. It also means that Fourier series expansion
technique could be applied on this application field.
A
numerical procedure to find the coefficients of Fourier series from a
series of given data was implemented in ABC FORTH. I did this
researching. Then I understood how to make a computer generate
different music sound. After many different kinds of sound were
generated via the specified code from the information manual, I found
that synthesized waveform was the main foundation. Other programming
factors, for example specifying lengths in time, mixing mode,
amplitude, frequency…etc, will mere influence effect of sound
output but will not influence the sound style. Piano was still piano,
guitar was still guitar, if waveform did not varied. Artificial
speech outputs are in the same condition.
To
simply the explanation for such an application, this article selected
a simple triangle waveform to prepare its data set. Directed Fourier
series formulas were used in the program. The coefficients of Fourier
series will be calculated out automatically. Data generator program
of outcome included as well for going to check. A practical
comparative diagram showed the result.
1. Question
Any
kind of periodic signal waveform could be represented by its
equivalent Fourier series, its general expression are:
f(x)
= ½ a(0) + a(1) cos(x) + a(2) cos(2x) + … + a(n) cos(nx)
+
b(1) sin(x) + b(2) sin(2x) + … + b(n) sin(nx) (1)
a(n)
= 1/pi ∫ f(x) cos(nx) dx
(2)
b(n)
= 1/pi ∫ f(x) sin(nx)
dx (3)
Integration
range is from 0 to 2pi.
(2), (3) expressions could be rewritten in terms of summation format:
a(n)
= 1/m ∑ f(k) cos ( nk * (2pi / 2m) ) (4)
b(n)
= 1/m ∑ f(k) sin ( nk * (2pi / 2m) ) (5)
Summation
range is from k = 1 to k = 2m.
According
to expressions (1) to (5), assume that k sets data were given, and m
= k/2. We are going to find n pairs of coefficients: a(n) and b(n), n
= k/2 as well.
Owing
to the integration in one cycle = 0, we got a(0) = 0, a(0) will be
neglected.
a(0)
= 1/pi ∫ f(x) dx = 0
Integration
range is from 0 to 2pi.
To
simplified our question, we restricted given data are exactly in 2pi periodic range only, and let 12 data points are distributed in equal
distance. Data are:
2pi*(1/12)
0.6666
2pi*(2/12)
0.3333
2pi*(3/12)
0.0
2pi*(4/12)
-0.3333
2pi*(5/12)
-0.6666
2pi*(6/12)
-1.0
2pi*(7/12)
-0.6666
2pi*(8/12)
-0.3333
2pi*(9/12)
0.0
2pi*(10/12)
0.3333
2pi*(11/12)
0.6666
2pi*(12/12)
1.0
The
waveform of given data is showing in fig.1.
Fig.1
What
are the coefficients of a(1) to a(6), b(1) to b(6)?
2. Code
\
Calculating Fourier coefficients in Lina64 ABC forth
\
===============================================
\
(1) Data structure declare heap
12
value km
km
2 / value nm
2
integers k n
km
array f
km
array t
km
array u
nm
array a
nm
array b
5
reals c d q r m
\
===============================================
\
===============================================
\
(2) Data input heap
:
input
{{
m = i>r ( nm ) }}
{{
f ( 0 ) = 1.0 e 0 }}
{{
f ( 1 ) = 6.666 e -1 }}
{{
f ( 2 ) = 3.333 e -1 }}
{{
f ( 3 ) = 0.0 e 0 }}
{{
f ( 4 ) = -3.333 e -1 }}
{{
f ( 5 ) = -6.666 e -1 }}
{{
f ( 6 ) = -1.0 e 0 }}
{{
f ( 7 ) = -6.666 e -1 }}
{{
f ( 8 ) = -3.333 e -1 }}
{{
f ( 9 ) = 0.0 e 0 }}
{{
f ( 10 ) = 3.333 e -1 }}
{{
f ( 11 ) = 6.666 e -1 }}
{{
f ( 12 ) = 1.0 e 0 }}
;
\
All parts above are changeable setting.
\
===============================================
\
===============================================
\
All parts below are fixed.
\
(3) Main code heap
:
main
basic
10
run input
20
for n = 1 to nm
30
let { q = 0.0 e 0 }
40
let { r = 0.0 e 0 }
50
for k = 1 to km
60
let { c = i>r ( k ) }
70
let { d = i>r ( n ) }
80
let { t ( k ) = ( f ( k ) * cos ( ( d * c * f(pi) ) / m ) ) / m }
90
let { u ( k ) = ( f ( k ) * sin ( ( d * c * f(pi) ) / m ) ) / m }
100
let { q = q + t ( k ) }
110
let { r = r + u ( k ) }
120
next k
130
let { a ( n ) = q }
140
let { b ( n ) = r }
150
print " format 1================== "
160
print " ( " ; n ; " ):: " , { q , r }
200
print " format 2=================== "
210
print " a ( " ; n ; " ) = " ; { q }
220
print " b ( " ; n ; " ) = " ; { r }
300
print " format 3=================== "
310
run cr ." a ( " n . ." ) = " q f.
320
run cr ." b ( " n . ." ) = " r f. cr cr
400
next n
410
end ;
\
===============================================
\
===============================================
\
(4) Manipulation heap
:
check
main
basic
20
for k = 0 to 360
30
let { c = i>r ( k ) }
40
let { d = f(2pi) * c / 360.0 e 0 }
50
let { m = a ( 1 ) * cos ( d ) + a ( 3 ) * cos ( 3.0 e 0 * d ) + a ( 5
) * cos ( 5.0 e 0 * d ) }
60
run cr d f. 5 spaces m f.
100
next k
110
end ;
\
================================================
3. Outcomes
main
format
1==================
(
1 ):: 8.2929502277 E -1 -2.484 E -17
format
2===================
a
( 1 ) = 8.2929502277 E -1
b
( 1 ) = -2.484 E -17
format
3===================
a
( 1 ) = 0.8292950227
b
( 1 ) = -0.0
format
1==================
(
2 ):: -3.0 E -19 -1.1 E -18
format
2===================
a
( 2 ) = -3.0 E -19
b
( 2 ) = -1.1 E -18
format
3===================
a
( 2 ) = -0.0
b
( 2 ) = -0.0
format
1==================
(
3 ):: 1.1113333333 E -1 -2.3 E -18
format
2===================
a
( 3 ) = 1.1113333333 E -1
b
( 3 ) = -2.3 E -18
format
3===================
a
( 3 ) = 0.1111333333
b
( 3 ) = -0.0
format
1==================
(
4 ):: -1.0 E -19 -2.0 E -19
format
2===================
a
( 4 ) = -1.0 E -19
b
( 4 ) = -2.0 E -19
format
3===================
a
( 4 ) = -0.0
b
( 4 ) = -0.0
format
1==================
(
5 ):: 5.9571643891 E -2 7.37 E -18
format
2===================
a
( 5 ) = 5.9571643891 E -2
b
( 5 ) = 7.37 E -18
format
3===================
a
( 5 ) = 0.0595716438
b
( 5 ) = 0.0
format
1==================
(
6 ):: 0.0 E 0 0.0 E 0
format
2===================
a
( 6 ) = 0.0 E 0
b
( 6 ) = 0.0 E 0
format
3===================
a
( 6 ) = 0.0
b
( 6 ) = 0.0
OK
4. Discussion
Despite
to calculate coefficients of Fourier series by numerical method is
not a complicated problem, owing to the calculation must cover on
full domain of sine and cosine, this question is a very nice
evaluation topic for testing of sine and cosine function. I
understood this merit. Before I consolidated my program, I had
modified many places in my ABC FORTH system to fit the application of
this question.
This
topic is totally different from the Fast Fourier Transform topic. All
analysis are still kept in the time domain, we don’t
have to consider frequency domain problem. Tedious calculation makes
us to do programming. In fact, this example question is an analytical
problem. You are able to find the analytical solution from a general
Fourier analysis textbook. Analytical solving procedures are tedious
as well. Numerical solution is relatively a directed and a simple
calculation.
In
this program, only one value km is a must to adjust variable. Km must
to be an even number. You are able to set out km = 100, then input
100 data, to get 50 pairs of a(n) and b(n). In real world,
engineering data are seldom could be represented as a simple
mathematical function, especially in speech record signal aspects.
This program is dealing with such kinds of input data.
Input
data in the question had been arranged in one period of cycle from 0
to 2pi,
and had been arranged symmetric to y-axis. Both conditions are
considered for simplify coding.
For
the sake of linking the using between my FORTH system and Python
Plotting system service, I had spent many days redesigned my floating
point number output instructions “F.”
and “FS.”.
At the beginning, both instructions will output total digits of
floating point number. After this article has been consolidated, both
instructions will output 12 digits (include “.”
or “-”
marks) at deferred setting condition. Variable setting of SIGFP can
help to adjust output digits. Besides, I adopted ignored treatment
for an infinitesimal number, if the number is less than 1.0E-18, or
if the number is less than the current setting 12 digits, F. output
will be 0.0.
In
Python plotting system, a number with decimal point but without E
mark inside will be considered as a floating point. This convenience
let me easy to design F. instruction in FORTH system. Now, all output
data are able to transfer to Python plotting program directly,
whatever how many sets. The following program generated 360 pairs of
data. According to it, Python plotting program will plot out the
overlapping diagram on the original triangle diagram. Fig.2 is
showing the achievement of my program.
:
check
main
basic
20
for k = 0 to 360
30
let { c = i>r ( k ) }
40
let { d = f(2pi) * c / 360.0 e 0 }
50
let { m = a ( 1 ) * cos ( d ) + a ( 3 ) * cos ( 3.0 e 0 * d ) + a ( 5
) * cos ( 5.0 e 0 * d ) }
60
run cr d f. 5 spaces m f.
100
next k
110
end ;
Fig.2
5. References
Hwei
P. Hsu. Ph. D. “Outline
of Fourier Analysis”,
Electrical Engineering department, Wayne State University, 1972
Chih-Fan
Chen, I. John Haas, “Elements
of Control Systems Analysis: Classical and Modern approaches”,
University of Houston, 1972
Erwin
Kreyszig, “Advanced
Engineering Mathematics”,
John Wiley and Sons, Inc. New York, 1967
Murray
R. Spiegel, “Fourier
Analysis”,
McGraw-Hill Book Company, 1974
6. My personal testing remarks
Instruction “F.” and “FS.”
approved.
沒有留言:
張貼留言