result of FFT

[Timon Gehr]

On 7/8/25 20:11, Matthew wrote:
> Hi,
> 
> I'm writing a program where I'm trying to decode DTMF tones. I already 
> completed the wave file decoder and I know I'm supposed to use an FFT to 
> transform the samples from time domain to frequency domain but I'm stuck 
> on determining which of the DTMF frequencies are present.
> 
> Consider the following, the input is the sound wave samples as an array 
> of floats between -1 and +1.
> 
> ```D
> void decode_sound(float[] samples)
> {
>      import std.numeric;
>      import std.math;
>      import std.complex;
>      import std.algorithm;
> 
>      size_t fft_size = 4096;
> 
>      auto f = new Fft(fft_size);
> 
>      foreach(ch; samples.chunks(fft_size))
>      {
>          auto res = f.fft(ch);
>          // res is an array of 4096 complex numbers
>      }
> }
> ```
> 
> I can't figure out how the 4096 results of the FFT relate to the 
> frequencies in the input.
> 
> I tried taking the magnitude of each element, I tried taking just the 
> real or imaginary parts.  I plotted them but the graph doesn't look how 
> I'm expecting.
> 
> What do the 4096 resulting complex numbers represent?
> How should I use the result to check whether the 1209Hz, 1336Hz, 1477Hz, 
> or 1633Hz tones are present in that part of the sound?
> 
> Thanks,
> Matthew

Well, with `N=fft_size`, `fft` computes

```
output[j]=∑ₖe^(-2πijk/N) input[k].
```

and we know that ifft gives:

```
input[k]=∑ⱼe^(2πijk/N) output[j]/N.
```

You are trying to write your input as a linear combination of sine 
waves, so we interpret `k` as time and `j/N` as frequency:

```
input[k]=∑ⱼe^(2πi k j/N)·output[j]/N.
```

Or, with `f_j[k]=e^(2πi k j/N)`:

```
input[k]=∑ⱼ(output[j]/N)·f_j[k].
```

One complication here is that both `f_j` and `output/N` are complex 
functions.

Note the relationship:

```
cos(ω·t+φ) = ½·(e^(i(ωt+φ))+e^(-i(ωt+φ))) = ½·(e^(iωt)·e^φ+e^(-iωt)·e^-φ)
```

In our setup, e^(iωt) corresponds to f_j, e^(-iωt) corresponds to 
f_(N-j), φ corresponds to the phase component of our outputs.

Interestingly, as your `input` signal is real, we actually know that 
`output[j]` and `output[N-j]` must be conjugate to each other as this is 
the only way to cancel out the imaginary parts.

This also implies that `output[N/2]` is real, which works out because 
`f_(N/2)[k]=(-1)^k`.


Putting it in practice:

```d
import std;

alias K=double;

enum sample_rate=44100;

auto fToWave(R)(size_t N,R coefficients_f){
     return iota(N/2+1).map!(j=>
         tuple!("magnitude","frequency","phase")(
             (j==N/2?1.0:2.0)*abs(coefficients_f[j]).re,
             K(j)*sample_rate/N,
             std.complex.log(coefficients_f[j]).im)
         );
}

void main(){
     enum N=4096;
     auto input=wave(N,1336.0f);
     auto output=fft(input);
     auto coefficients_f=output.map!(o=>o/N);
     //auto input2=linearCombination(N,coefficients_f,f(N)); // (slow, 
just for illustration)
     //assert(isClose(input,input2,1e-7,1e-7));
     auto coefficients_w=fToWave(N,coefficients_f);
     //auto input3=reconstructWave(N,coefficients_w); // (slow, just for 
illustration)
     //assert(isClose(input,input3,1e-3,1e-3));
     auto sorted=coefficients_w.array.sort.retro;
     foreach(magnitude,frequency,phase;sorted.until!(coeff=>coeff[0]<1e-2)){
  
writefln!"%s·cos(2π·%s·t%s%s)"(magnitude,frequency,text(phase).startsWith("-")?"":"+",phase);
     }
}

// helper functions for illustration

auto time(size_t N)=>iota(N).map!(i=>K(i)/sample_rate);

auto wave(size_t N,K frequency,K phase=0.0f){
     K angularVelocity=2*PI*frequency;
     return time(N).map!(t=>cos(angularVelocity*t+phase));
}

/+
// helper functions for slow illustrations:

enum I=Complex!K(0,1);

auto f(size_t N){
     return iota(N).map!(j=>iota(N).map!(k=>exp(2*PI*I*j*k/N)));
}

auto linearCombination(R,RoR)(size_t N,R coefficients,RoR functions){ // 
compute ∑ⱼcoefficients[j]·functions[j]
     return 
iota(N).map!(k=>sum(iota(coefficients.length).map!(j=>coefficients[j]*functions[j][k])));
}

auto reconstructWave(R)(size_t N,R coefficients_w){
     auto coefficients=coefficients_w.map!(coeff=>coeff.magnitude);
     auto 
functions=coefficients_w.map!(coeff=>wave(N,coeff.frequency,coeff.phase));
     return linearCombination(N,coefficients,functions);
}
+/

```

(Note that additional processing may be helpful, as `fft` is based on 
the assumption that signals repeat with a period of `N`.)