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`.)
