Fast Fourier transform: Difference between revisions

From APL Wiki
Jump to navigation Jump to search
(→‎Simple recursive implementation: otherwise produces text garbage on incorrect input)
m (Text replacement - "<source" to "<syntaxhighlight")
 
(One intermediate revision by the same user not shown)
Line 5: Line 5:
== Implementations ==
== Implementations ==
=== APLX ===
=== APLX ===
This FFT code is implemented with the [[wikipedia:Cooley–Tukey FFT algorithm|Cooley–Tukey FFT algorithm]] by dividing the transform into two pieces of size <source lang=apl inline>N÷2</source> at each step. It will run under [[APLX]].
This FFT code is implemented with the [[wikipedia:Cooley–Tukey FFT algorithm|Cooley–Tukey FFT algorithm]] by dividing the transform into two pieces of size <syntaxhighlight lang=apl inline>N÷2</syntaxhighlight> at each step. It will run under [[APLX]].


This is as given in Robert J. Korsan's article in APL Congress 1973, p 259-268, with just line labels and a few comments added.
This is as given in Robert J. Korsan's article in APL Congress 1973, p 259-268, with just line labels and a few comments added.


* X and Z are two-row [[Matrix|matrices]] representing the input and output real and imaginary data. The data length must be <source lang=apl inline>2*N</source> (N integer), and the algorithm will cope with varying N, unlike many DSP versions which are for fixed N.
* X and Z are two-row [[Matrix|matrices]] representing the input and output real and imaginary data. The data length must be <syntaxhighlight lang=apl inline>2*N</syntaxhighlight> (N integer), and the algorithm will cope with varying N, unlike many DSP versions which are for fixed N.
* Zero frequency is at <source lang=apl inline>Z[1;]</source>, maximum frequency in the middle; from there to <source lang=apl inline>¯1↑[1] Z</source> are negative frequencies. i.e. for an input Gaussian they transform a 'bath-tub' to a 'bath-tub'.
* Zero frequency is at <syntaxhighlight lang=apl inline>Z[1;]</syntaxhighlight>, maximum frequency in the middle; from there to <syntaxhighlight lang=apl inline>¯1↑[1] Z</syntaxhighlight> are negative frequencies. i.e. for an input Gaussian they transform a 'bath-tub' to a 'bath-tub'.
* This is an elegant algorithm, and works by transforming the input data into an array of 2×2 [[wikipedia:Butterfly diagram|FFT Butterflies]].
* This is an elegant algorithm, and works by transforming the input data into an array of 2×2 [[wikipedia:Butterfly diagram|FFT Butterflies]].


<source lang=apl>
<syntaxhighlight lang=apl>
     Z←fft X;N;R;M;L;P;Q;S;T;O
     Z←fft X;N;R;M;L;P;Q;S;T;O
⍝
⍝
Line 30: Line 30:
→loop
→loop
done:Z←⍉(N,2)⍴(+⌿X),[O-0.5]-⌿X
done:Z←⍉(N,2)⍴(+⌿X),[O-0.5]-⌿X
</source>
</syntaxhighlight>


=== Simple recursive implementation ===
=== Simple recursive implementation ===
<source lang=apl>
<syntaxhighlight lang=apl>
fft←{
fft←{
     1>K←2÷⍨M←≢⍵:⍵
     1>K←2÷⍨M←≢⍵:⍵
Line 42: Line 42:
     (odd+T),odd-T←even×exp
     (odd+T),odd-T←even×exp
}
}
</source>
</syntaxhighlight>
{{Works in|[[Dyalog APL]]}}
{{Works in|[[Dyalog APL]]}}
Sample usage:
Sample usage:
Line 50: Line 50:


Inverse FFT can be defined for testing:
Inverse FFT can be defined for testing:
<source lang=apl>
<syntaxhighlight lang=apl>
       ifft←{(≢⍵)÷⍨+fft+⍵}
       ifft←{(≢⍵)÷⍨+fft+⍵}
       test←{⌈/(10○⊢)(⍵-ifft fft ⍵)}
       test←{⌈/(10○⊢)(⍵-ifft fft ⍵)}
       test 1 1 1 1 0 0 0 0
       test 1 1 1 1 0 0 0 0
7.850462E¯17
7.850462E¯17
</source>
</syntaxhighlight>


2-dimensional FFT and inverse 2D FFT:
2-dimensional FFT and inverse 2D FFT:
<source lang=apl>
<syntaxhighlight lang=apl>
fft2D←{
fft2D←{
     ∨/0≠1|2⍟⍴⍵:'Matrix dimensions must be powers of 2'
     ∨/0≠1|2⍟⍴⍵:'Matrix dimensions must be powers of 2'
Line 64: Line 64:
}
}
ifft2D←{(≢∊⍵)÷⍨+fft2D+⍵}
ifft2D←{(≢∊⍵)÷⍨+fft2D+⍵}
</source>
</syntaxhighlight>


Sample usage:
Sample usage:
Line 76: Line 76:
=== Dyalog APL ===
=== Dyalog APL ===
FFT appears in [[dfns.dws]], a [[workspace]] supplied with [[Dyalog APL]], in the context of fast multi-digit multiplication<ref>dfns.dws: [http://dfns.dyalog.com/n_xtimes.htm xtimes] — Fast multi-digit product using FFT</ref>. Extracted from there, it is there defined as:
FFT appears in [[dfns.dws]], a [[workspace]] supplied with [[Dyalog APL]], in the context of fast multi-digit multiplication<ref>dfns.dws: [http://dfns.dyalog.com/n_xtimes.htm xtimes] — Fast multi-digit product using FFT</ref>. Extracted from there, it is there defined as:
<source lang=apl>
<syntaxhighlight lang=apl>
roots←{×\1,1↓(⍵÷2)⍴¯1*2÷⍵}
roots←{×\1,1↓(⍵÷2)⍴¯1*2÷⍵}
cube←{⍵⍴⍨2⍴⍨2⍟⍴⍵}
cube←{⍵⍴⍨2⍴⍨2⍟⍴⍵}
floop←{(⊣/⍺)∇⍣(×m)⊢(+⌿⍵),[m-0.5]⍺×[⍳m←≢⍴⍺]-⌿⍵}
floop←{(⊣/⍺)∇⍣(×m)⊢(+⌿⍵),[m-0.5]⍺×[⍳m←≢⍴⍺]-⌿⍵}
FFT←{,(cube roots⍴⍵)floop cube ⍵}
FFT←{,(cube roots⍴⍵)floop cube ⍵}
</source>
</syntaxhighlight>


== References ==
== References ==

Latest revision as of 21:06, 10 September 2022

The fast Fourier transform (FFT) is an algorithm to compute the discrete Fourier transform of a vector in time , where a naive implementation achieves only time. APL implementations of the fast Fourier transform began appearing as early as 1970, with an 8-line implementation by Alan R. Jones published in APL Quote-Quad.[1]

A Fourier Transform (FFT) is a method of calculating the frequency components in a data set — and the inverse FFT converts back from the frequency domain — 4 applications of the FFT rotates you round the complex plane and leaves you back with the original data.

Implementations

APLX

This FFT code is implemented with the Cooley–Tukey FFT algorithm by dividing the transform into two pieces of size N÷2 at each step. It will run under APLX.

This is as given in Robert J. Korsan's article in APL Congress 1973, p 259-268, with just line labels and a few comments added.

  • X and Z are two-row matrices representing the input and output real and imaginary data. The data length must be 2*N (N integer), and the algorithm will cope with varying N, unlike many DSP versions which are for fixed N.
  • Zero frequency is at Z[1;], maximum frequency in the middle; from there to ¯1↑[1] Z are negative frequencies. i.e. for an input Gaussian they transform a 'bath-tub' to a 'bath-tub'.
  • This is an elegant algorithm, and works by transforming the input data into an array of 2×2 FFT Butterflies.
    Z←fft X;N;R;M;L;P;Q;S;T;O
⍝
⍝ Apl Congress 1973, p 267. Robert J. Korsan.
⍝
⍝ Restructure as an array of primitive 2×2 FFT Butterflies
X←(2,R←(M←⌊2⍟N←¯1↑⍴X)⍴2)⍴⍉X
⍝ Build sin and cosine table :
Z←R⍴⍉2 1∘.○○(-(O←?1)-⍳P)÷P←N÷2
⍝
Q←⍳P←M-1+L←0
T←M-~O
loop:→(M≤L←L+1)⍴done
X←(+⌿X),[O+¯0.5+S←M-L](-/Z×-⌿X),[O+P-0.5]+/Z×⌽-⌿X
Z←(((-L)⌽Q),T)⍉R⍴((1+P↑(S-1)⍴1),2)↑Z
→loop
done:Z←⍉(N,2)⍴(+⌿X),[O-0.5]-⌿X

Simple recursive implementation

fft←{
    1>K←2÷⍨M←≢⍵:⍵
    0≠1|2⍟M:'Length of the input vector must be a power of 2'
    odd←∇(M⍴1 0)/⍵
    even←∇(M⍴0 1)/⍵
    exp←*(0J¯2×(○1)×(¯1+⍳K)÷M)
    (odd+T),odd-T←even×exp
}
Works in: Dyalog APL

Sample usage: 

      fft 1 1 1 1 0 0 0 0
4 1J¯2.414213562 0 1J¯0.4142135624 0 1J0.4142135624 0 1J2.414213562

Inverse FFT can be defined for testing: 

      ifft←{(≢⍵)÷⍨+fft+⍵}
      test←{⌈/(10○⊢)(⍵-ifft fft ⍵)}
      test 1 1 1 1 0 0 0 0
7.850462E¯17

2-dimensional FFT and inverse 2D FFT: 

fft2D←{
    ∨/0≠1|2⍟⍴⍵:'Matrix dimensions must be powers of 2'
    ⍉(fft⍤1)⍉(fft⍤1)⍵
}
ifft2D←{(≢∊⍵)÷⍨+fft2D+⍵}

Sample usage: 

      fft2D 2 2⍴⍳4
10 ¯2
¯4  0
      ifft2D fft2D 2 2⍴⍳4
1 2
3 4

Dyalog APL

FFT appears in dfns.dws, a workspace supplied with Dyalog APL, in the context of fast multi-digit multiplication[2]. Extracted from there, it is there defined as: 

roots←{×\1,1↓(⍵÷2)⍴¯1*2÷⍵}
cube←{⍵⍴⍨2⍴⍨2⍟⍴⍵}
floop←{(⊣/⍺)∇⍣(×m)⊢(+⌿⍵),[m-0.5]⍺×[⍳m←≢⍴⍺]-⌿⍵}
FFT←{,(cube roots⍴⍵)floop cube ⍵}

References

  1. ↑ Jones, Alan R. (IBM). "Fast Fourier transform". APL Quote-Quad Volume 1, Number 4. 1970-01.
  2. ↑ dfns.dws: xtimes — Fast multi-digit product using FFT
Retrieved from ‘https://aplwiki.com/index.php?title=Fast_Fourier_transform&oldid=9195’