Timon Gehr 
Posted in reply to Matthew
| 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`.)
|
Permalink
Reply
