This is an relatively unknown Fourier method that is very suitable
for generating harmonic amplitudes on a sample-by-sample basis.
As such it is interesting for use in audio tools and circuit
simulation. The method can compete with the FFT, especially when
only a few harmonics are needed and the base frequency is known.

I have to admit that I initially coded this in MATLAB. However,
only when I redid it again in Forth the details became really
clear. The iForth code is 5 times faster than MATLAB.

The results are visually checked through Gnuplot. iForth provides
a MATLAB-like interface to that program. I've added a link to
a pdf on my Github repository.

-marcel

-- --------------------------------------------------
https://github.com/mhx-gh/SYSSIM/blob/master/mft.pdf
-- --------------------------------------------------
(*
* LANGUAGE : ANS Forth with extensions
* PROJECT : Forth Environments
* DESCRIPTION : The Momentary Fourier Transform
* CATEGORY : Tool for circuit simulation
* AUTHOR : Marcel Hendrix
* LAST CHANGE : 21:05:01, January 22, 2024, Marcel Hendrix
*)

NEEDS -miscutil
NEEDS -gaussj
NESTING @ 1 =
[IF] NEEDS -gnuplot
[THEN]
REVISION -MFT "--- Momentary FFT Version 0.03 ---"

PRIVATES

DOC
(*
S. Albrecht, I. Cumming and J. Dudas, "The momentary Fourier
transformation derived from recursive matrix transformations,"
Proceedings of 13th International Conference on Digital Signal
Processing, Santorini, Greece, 1997, pp. 337-340 Vol. 1.
Compute a maximum of N harmonics from a signal's N-point sample
record. Complex arithmetic is needed.

The nice thing about the algorithm is that, because the result is
build sample by sample, one doesn't need to have the (equidistant)
samples in a block of size N.
When sending a repeating waveform with N samples per cycle, a plot of
the result shows the amplitude of the N harmonics changing during the
first cycle, resulting in the expected values at the end of that cycle.
From that point onwards the amplitude of the harmonics stays constant.
The plot shows transients when the signal changes between blocks
of N samples.

For use in iSPICE, preface with a SAMPLEr (not a samplehold!) that
triggers (MFT) for every sample taken. The output of the MFT subcircuit
can contain upto N output nodes providing the value of the harmonic
amplitudes. The *phases* of these harmonics can be displayed, however,
they change N times per basic interval (by 360/N degrees). A way to
decide on a meaningful start of the sample block is needed when
requiring phase.
*)
ENDDOC

0 VALUE N PRIVATE -- how many samples there are in a cycle
DOUBLE DARRAY x{ PRIVATE -- the input sample buffer
COMPLEX DARRAY y{ PRIVATE -- contains N Fourier coefficients of x{.
COMPLEX DARRAY twiddle{ PRIVATE -- turns the transformation into a DFT

-- The FFT twiddle factors are: twiddle(k) = exp(-2*pi*k/N)
-- but here we want 1 / twiddle(k) = exp(2*pi*k/N).
-- The paper makes this a matrix, but it is faster as a vector.
: (MFT) ( N -- y{ ) ( F: xin -- )
TO N
y{ CDIM N \ we want to re-initialize when N changes
<> IF x{ N }malloc
y{ N }malloc
twiddle{ N }malloc malloc-fail? ?ALLOCATE
( xin ) x{ N 1- } DF!
PI*2 N S>F F/
N 0 ?DO FDUP I S>F F* FEXP(jX) twiddle{ I } Z! LOOP
FDROP
y{ EXIT
ENDIF
\ apply the MFT to sample xin
x{ 0 } DF@ FLOCAL sos \ shifted-out sample
x{ 1 } x{ 0 } N 1- DFLOATS MOVE \ circshift (x, LEFT)
( xin ) x{ N 1- } DF!
N 0 ?DO
x{ N 1- } DF@
sos F- REAL>Z \ (xin - prev + y(i)) * twiddle(i)
y{ I } DUP Z@ Z+
twiddle{ I } Z@ Z*
( addr -- ) Z!
LOOP y{ ; \ needs scaling by 2/N & DC term/2

NESTING @ 1 =
[IF] ( plotting with Gnuplot )

100e FVALUE fin
#64 VALUE Ns
2 VALUE m

m Ns * DOUBLE ARRAY t{
m Ns * DOUBLE ARRAY s{
Ns m Ns * DOUBLE MATRIX f{{

fin Ns S>F F* FVALUE fs
fs 1/F FVALUE Ts

: INIT ( -- ) \ construct a signal with H0=0.1, H1=0.8, H2=0.6, H3=0.4
0e FLOCAL time
m Ns * 0 ?DO Ts I S>F F* FDUP TO time t{ I } DF!
0.1e
fin PI*2 F* time F* FSIN 0.8e F* F+
fin F2* PI*2 F* time F* FSIN 0.6e F* F+
fin 3e F* PI*2 F* time F* FSIN 0.4e F* F+
s{ I } DF!
LOOP ;

: GET-RESULTS ( -- )
2e Ns S>F F/ FLOCAL 2/N
m Ns * 0 ?DO s{ I } DF@ Ns (MFT) ( y{ -- )
>S Ns 0 ?DO S I } Z@ ZABS
2/N F*
I 0= IF F2/ ENDIF ( dc )
f{{ I J }} DF!
LOOP -S
LOOP ;

: PLOT-RESULTS ( -- )
0 FIGURE CLF
S" MFT shell test program" TITLE
2 1 SUBPLOT
S" amplitude [V]" YLABEL
TRUE GRID TRUE ZOOM
[PHOLD t{ s{ S" blue-" PLOT PHOLD]
S" amplitude [%]" YLABEL
S" time [ms]" XLABEL
TRUE GRID TRUE ZOOM
[PHOLD
0 TO subp
S" dc" S" h1" S" h2" S" h3" -4 LEGEND
t{ f{{ 0 =rowget 100e =scalemat S" black-" PLOT
t{ f{{ 1 =rowget 100e =scalemat S" red-" PLOT
t{ f{{ 2 =rowget 100e =scalemat S" blue-" PLOT
t{ f{{ 3 =rowget 100e =scalemat S" magenta-" PLOT
PHOLD]
SPEND ;

: MFT_shell ( -- ) INIT GET-RESULTS PLOT-RESULTS ;

:ABOUT CR ." Try: MFT_shell"
CR ." Expect DC = 10%, H1=80%, H2=60%, H3=40% " ;

.ABOUT -MFT
[THEN]

DEPRIVE

(* End of Source *)