by Theo Verelst
Using the Tcl/Tk package Bwise, and some fundamental, tcl-coded, signal processing procedures, such as fourier transformation, I'll examplify the important signal processing concept of frequency modulation by making an Tk applet to play around with this modulation principle, known in radio engineering, wireless communication, modems, and sound synthesis, to mention a few fields.
The interest comes mainly from two signal processing considerations (are the signal processing page leaders on holiday ?! ;) ), first, it is an overlooked principle/theorym which is of great value in many signal processing/analysis applications, also when tcl is used (as appears regulalrly on these pages), and second, I think it is a fun and interesting and relevant use of tcl/tk to examplify an important principle, which for me historically goes back to the advent of the Yamaha DX-7 synthesizer which is based on sine wave frequency modulation sound synthesis. When I was a second year electrical engineering student looking for a section to graduate in, I had just worked enough to buy a DX-7, and I had the schematic diagram of it, because I was very interested in what the chips (at the time competatively big and complicated ones) in that synthesizer worked like, and also how FM synthesis worked, so at the time I programmed a sample reproducer and FM program on the Atari ST. At the time, such programs were preferably compiled, but nowadays we can use Tcl and do things more or less in interaction time, and even make an interactive user interface in Tk on top. IMNSHO a very important, and protection-worthy use of a modern progamming language like tcl/tk.
First, lets make a few one-liners to de signal frequency analysis, first a 256 (tcl-list) sample sinewave (three wavelengths are fitted in the list in this case, so a third harmonic):
set twopi 2.0*3.1415926535
set r {}; for {set i 0} {$i < 256} {incr i} {lappend r [expr sin( 3.0 * $twopi*$i/256)] }To find a certain harmonic in the signal, we can use the well known fourier analysis, in this case in discrete form (multiplying the signal sample for sample with a certain harmonic component):
set t 0.0; for {set i 0} {$i < 256} {incr i} {set t [expr $t+ [lindex $r $i] * sin( 3.0 * $twopi*$i/256)] } ; puts [ expr $t/([llength $r] /2.0)]The third harmonic returns a 'strength' of 1 (full amplitude), then the 3.0 above is replaced by 1.0 for fundamantal (1st harmonic) or another harmonic (smaller than 128) some number very close to zero (background numerical noise) is returned.
I'd rather make an application which does something like this [L1 ] (Click on the left on Waveform Lab, Theory, View Applet, Making Wavews, Sine f=2) in tcl/tk, but I want to do other things first, and I don't know if something like it already exists.
Lets first make sure we can get some of the important spectral lines from a tcl-list of samples, first, a fourier transform proc based on the above oneliner tries, and a filter for all components higher then 1e-10:
proc fourier { {s} } {
set r {}
set twopi 2.0*3.1415926535
for {set j 0} {$j < [expr [llength $s]/2]} {incr j} {
set t 0.0
for {set i 0} {$i < 256} {incr i} {set t [expr $t+ [lindex $s $i] * sin( $j * $twopi*$i/256)] };
lappend r [ expr $t/([llength $s] /2.0)]
}
return $r
}
set i 0; foreach j [fourier $r] {if {[expr abs($j)] > 1e-10} {puts "$i $j"} ; incr i}I've used the signal from list r, which can contain any number of samples with harmonics of a multiple of largest wave which fits in the interval, up to half te number of samples -1. Note that I've not used an FFT (Fast Fourier Transform, but just a convolution of the signal wil each of the possible harmmonics in the interval seperately.
Lets try a signal with more components before we plot it and the discrete fourier transform result in a Tk canvas, and compare the result with the target: FM modulates signals, which we'll also store in tcl lists, as sampled signals. Note that The above fourier routine only takes what would be called the imaginary component, only the sine component, which in fact we should take both the sine and cosine components into account, also when we compute the norm of each component.
(TV aug 10, 2004, 19:47 Ooops, (local) power outage in Amsterdam, I'll see what I can do right now, might be I have to get some car power as soon as the notebook is running out of current... For the moment, mind that my images on the wiki are served from a not completely up to date server, though most is there. 20:01, ok all data should be served in latest version again!)
A few more harmonics:
set s {} ; for {set i 0} {$i < 256} {incr i} { set t 0.0; foreach {h a} {1 1 2 0.5 4 0.25} { set t [expr $t + $a*sin($h*$twopi*$i/256)]} ; lappend s $t }Testing:
set i 0; foreach j [fourier $s] {if {[expr abs($j)] > 1e-10} {puts "$i $j"} ; incr i}
1 1.00000000004
2 0.500000000013
4 0.250000000003(Beware this is not a fast fourier transform, you won't like it on a 75 MHz pentium, which runs tcl/tk just fine..)
Or using a procedure to add sine components together which it makes from harmonic, intensity tuples as argument:
proc additive { {c {1.0 1.0}} {l 256} } {
set twopi 2.0*3.1415926535
set s {} ;
for {set i 0} {$i < $l} {incr i} {
set t 0.0;
foreach {h a} $c {
set t [expr $t + $a*sin($h*$twopi*$i/$l)]
} ;
lappend s $t
}
return $s
}
set i 0; foreach j [fourier [additive {1 1 2 0.5 33 0.4411}]] {if {[expr abs($j)] > 1e-10} {puts "$i $j"} ; incr i}
1 1.00000000005
2 0.500000000013
33 0.441100000007(This reminds me of highschool where I'd make organ imitation waves like this in BASIC, which I would push into a DA converter with a machine language routine. Nowadays in principle snack could be used for the purpose.)
Ok, so that looks fine thus far, too.
We should be aware of our statistics, self correlation margin, or measurement accuracy here: I only take integer harmonics into account, which is the best one can do on a limited interval, with no repeated (and averaged) harmonics measurements, so for instance a harmonics at one third of maximum shannon frequency (1/6th of the sample frequency) gives clearly an unclean spectrum:
set i 0; foreach j [fourier [additive [list 1 1 2 0.5 33 0.4411 [expr 256.0 /6] [expr 1.0/3] ] ]] {if {[expr abs($j)] > 1e-10} {puts "$i $j"} ; incr i}
1 0.999944619477
2 0.4998890717
3 -0.000166812255595
4 -0.000223205146321 5 -0.000280285428205 6 -0.000338239357798 7 -0.000397263292805
8 -0.000457566242402 9 -0.000519372742715 10 -0.000582926134539 11 -0.000648492370484
12 -0.00071636445592 13 -0.000786867720594 14 -0.000860366104961 15 -0.000937269742203
16 -0.0010180441847 17 -0.00110322172359 18 -0.00119341541385 19 -0.00128933659812
20 -0.00139181702883 21 -0.00150183706734 22 -0.00162056203473 23 -0.0017493896104
24 -0.00189001243408 25 -0.00204450193346 26 -0.00221542229237 27 -0.00240598804451
28 -0.00262028610213 29 -0.00286359523334 30 -0.00314285683554 31 -0.00346738780974
32 -0.00384999438229 33 0.436791222761 34 -0.00487018638954 35 -0.00557446661084
36 -0.00648600400252 37 -0.00771462747441 38 -0.00946429792141 39 -0.0121614705416
40 -0.0168722722987 41 -0.0272215182033 42 -0.0685828738951 43 0.138153907998
44 0.0347683150063 45 0.0199892750175 46 0.0140709960594 47 0.0108791817661
48 0.00888025442656 49 0.00750920881886 50 0.00650922199082 51 0.00574674399131
52 0.00514545883687 53 0.00465859072406 54 0.00425587165083 55 0.00391684911001
56 0.00362720532238 57 0.00337661724629 58 0.00315745397813 59 0.0029639526553
60 0.00279167961981 61 0.00263716792751 62 0.00249766735467 63 0.00237096817087
64 0.00225527447239 65 0.00214911154935 66 0.00205125706634 67 0.00196068921301
68 0.00187654713692 69 0.00179810038757 70 0.0017247250688 71 0.00165588503881
72 0.00159111695732 73 0.00153001828748 74 0.00147223759734 75 0.00141746666083
76 0.00136543397471 77 0.00131589941348 78 0.00126864979138 79 0.00122349515452
80 0.00118026566841 81 0.0011388089946 82 0.00109898806346 83 0.00106067918278
84 0.0010237704173 85 0.000988160198281 86 0.000953756124445 87 0.000920473925813
88 0.00088823656525 89 0.000856973453391 90 0.000826619765383 91 0.000797115839992
92 0.000768406655256 93 0.000740441360005 94 0.000713172864626 95 0.000686557476632
96 0.000660554578524 97 0.000635126338873 98 0.000610237455044 99 0.000585854922419
100 0.000561947828895 101 0.000538487168927 102 0.000515445675413 103 0.000492797669845
104 0.000470518924218 105 0.00044858653717 106 0.00042697882192 107 0.000405675203593
108 0.000384656123691 109 0.000363902954213 110 0.000343397921173 111 0.000323124029927
112 0.000303064997844 113 0.000283205193617 114 0.000263529579362 115 0.000244023656798
116 0.000224673418637 117 0.000205465301751 118 0.000186386141168 119 0.000167423130528
120 0.000148563783388 121 0.000129795896384 122 0.000111107511816 123 9.24868858195e-005
124 7.3922454711e-005 125 5.54028034085e-005 126 3.69166343052e-005 127 1.84527380447e-005Ugly. If we'd take a multiple of this number samples, and assume the signal is repetetive, things would turn out better by averaging, but we'd need a lot of additional sample lists to come down to under 1e-10, it appears. Something to take into account when analyzing the spectrum of a bessel function (tcllib ??) gouverned Frequency Modulated signal. To check that we didn't just take our harmonic too high, lets take the 85th harmonic:
set i 0; foreach j [fourier [additive [list 1 1 2 0.5 33 0.4411 [expr 85.0 ] [expr 1.0/3] ] ]] {if {[expr abs($j)] > 1e-10} {puts "$i $j"} ; incr i}
1 1.00000000002
2 0.50000000001
33 0.441100000005
85 0.333333333309Ok.