Android 撥號音DTMF 編碼解碼理論和實戰
一、原理，公式推導
f(n)=2*cos(w*Ts)*f(n-1) - f(n-2)
float a1=2*cos(w*Ts);
float a2=-1;
假如我們需要產生取樣頻率為8KHz的440Hz的正弦波，那麼
a1=2*cos(2*pi*440/8000)=1.8817615，
而y[1]=sin(2*pi*440/8000)=0.33873792。
現在看如何用定點小數來更快的計算正弦波。我們使用16bit也就是short型的整數來表示定點小數。首先需要決定的是小數的Q值，雖然我們最後計算的正弦波的值都是小於1的，但是在計算過程中需要用2*cos(w*Ts)，而這個值最大為2，所以我們選擇的Q值必須至少最大能表示2。這裡我們選擇 Q14，Q14的定點小數
下面就是完整的計算8KHz取樣頻率的400Hz的定點小數的正弦波的程式。
以下內容為程式程式碼:
short y[3] = {0, 0x15AE,0}; // y(n), y(n-1), y(n-2)
short a1=0x786F;
short a2=0xC000;
short singen(){
y[0]=( (long)a1*(long)y[1]+(long)a2*(long)y[2] )>>14;
y[2]=y[1]; //// 遞迴推導公式，f(n)
y[1]=y[0];
return y[0];
}
使用定點小數計算不但速度比浮點更快，而且計算得出來的值是整數，這個數值可以直接傳遞給DAC(數模轉換器)轉換為模擬的聲音訊號，如果使用浮點小數計算的話，還必須把浮點數轉換為整數才能傳遞給DAC。
使用定點小數計算必須仔細分析誤差，下面來看看我們產生的正弦波的誤差是多少。定點小數計算中的誤差就是由定點小數表達精度決定的。在上面的例子中我們用 0x786F表示1.8817615，這存在一定的誤差，把Q14的0x786F再轉換為浮點數就是0x786F/2^14=1.8817749，我們可以看到相對誤差非常小，也就是說最終得到的正弦波在頻率上的誤差也是非常小的。
但是，定點小數並不是什麼時候都這麼精確。例如如果用CD 音質的取樣頻率44100Hz來產生100Hz的正弦波，那麼a1=2*cos(2*pi*440/44100)= 1.9960713，這個數轉換為16位元的Q14的值是0x7fc0。我們可以看到這時定點小數已經十分接近0x7fff了，最終產生的正弦波的頻率也會有很大的誤差。為了能夠精確地計算這樣的正弦波，必須使用32bit的Q30定點小數。關於32bit定點小數的計算方法將在別的章節介紹。
另外上面的singen函式每呼叫一次只產生一個值，如果要產生實時的正弦波的話，函式的呼叫頻率和取樣頻率相同，DSP的負擔相對比較大。一般DSP計算都採取塊計算方式，一次計算n個（例如64）個取樣值，這樣不但減少了函式的呼叫負擔，也可以減少中間的記憶體移動的次數（y[2]=y[1];y[1]= y[0];）。
上面的singen函式每呼叫一次只產生一個值，如果要產生實時的正弦波的話，函式的呼叫頻率和取樣頻率相同，DSP的負擔相對比較大。一般DSP計算都採取塊計算方式，一次計算n個（例如64）個取樣值，這樣不但減少了函式的呼叫負擔，也可以減少中間的記憶體移動的次數（y[2]=y[1];y[1]= y[0];）。
二、 對應Android程式碼分解
檔案： ./media/libaudioclient/ToneGenerator.cpp
const ToneGenerator::ToneDescriptor ToneGenerator::sToneDescriptors[] = {
{ .segments = {
{ .duration = ToneGenerator::TONEGEN_INF,
.waveFreq = { 1336, 941, 0 }, 0, 0},
{ .duration = 0 ,
.waveFreq = { 0 }, 0, 0}}, /////暫停時間，間隔
.repeatCnt = ToneGenerator::TONEGEN_INF,
.repeatSegment = 0 }, // TONE_DTMF_0
........................ // 看到這個就明白了
bool ToneGenerator::prepareWave() {
///// 準備好13個鍵的 Tone，f1， f2資料， 0 ~9 * #
while (mpToneDesc->segments[segmentIdx].duration) { //// 13 個 Segment
// Get total number of sine waves: needed to adapt sine wave gain.
unsigned int lNumWaves = numWaves(segmentIdx);
unsigned int freqIdx = 0;
unsigned int frequency = mpToneDesc->segments[segmentIdx].waveFreq[freqIdx];
while (frequency) {
// Instantiate a wave generator if ot already done for this frequency
if (mWaveGens.indexOfKey(frequency) == NAME_NOT_FOUND) {
ToneGenerator::WaveGenerator *lpWaveGen =
new ToneGenerator::WaveGenerator(mSamplingRate,
frequency,
TONEGEN_GAIN/lNumWaves);
mWaveGens.add(frequency, lpWaveGen);
//// 將頻率和物件存入向量中，以供後續使用
}
frequency = mpNewToneDesc->segments[segmentIdx].waveFreq[++freqIdx];
}
segmentIdx++;
void ToneGenerator::audioCallback(int event, void* user, void *info) {
////// 由於tone不是聲音檔案，需要按規則生成
unsigned int lFreqIdx = 0;
unsigned short lFrequency = lpToneDesc->segments[lpToneGen->mCurSegment].waveFreq[lFreqIdx];
while (lFrequency != 0) {
WaveGenerator *lpWaveGen = lpToneGen->mWaveGens.valueFor(lFrequency);
lpWaveGen->getSamples(lpOut, lGenSmp, lWaveCmd);
/// 請注意：lpOut，用的巧，高低頻率都在這裡加上了
lFrequency = lpToneDesc->segments[lpToneGen->mCurSegment].waveFreq[++lFreqIdx];
}
ToneGenerator::WaveGenerator::WaveGenerator(uint32_t samplingRate,
unsigned short frequency, float volume) {
F_div_Fs = frequency / (double)samplingRate;
d0 = - (float)GEN_AMP * sin(2 * M_PI * F_div_Fs);
mS2_0 = (short)d0;
mAmplitude_Q15 = (short)(32767. * 32767. * volume / GEN_AMP); // 這個不太理解
d0 = 32768.0 * cos(2 * M_PI * F_div_Fs); // Q15*2*cos()
void ToneGenerator::WaveGenerator::getSamples(short *outBuffer,
// loop generation
while (count) {
count--;
Sample = ((lA1 * lS1) >> S_Q14) - lS2;
// shift delay
lS2 = lS1;
lS1 = Sample;
Sample = (lAmplitude * Sample) >> S_Q15;
*(outBuffer++) += (short)Sample; // put result in buffer，使用加號很有意義
}
三、Dtmf音的識別，反解碼演算法
從某些連續段落的波形資料中找到是某兩個已定義的高頻與低頻的組合, 即可對應出 tone digit.
我們使用的是 goertzel 演算法,
這個演算法的效率比 DFT 或是 FFT 都高. 一組波形資料經過 Goertzel 演算法後都會得到該頻率的匹配度.
def goertzel_mag(numSamples, target_freq, sample_rate, data):
'''
int k, i
float floatnumSamples
float omega, sine, cosine, coeff, q0, q1, r2, magnitude, real, imag
'''
#floatnumSamples = (float)numSamples
scalingFactor = numSamples / 2.0
k = (int) (0.5 + ((numSamples * target_freq)/sample_rate))
omega = (2.0 * math.pi * k)/numSamples
sine = math.sin(omega)
cosine = math.cos(omega)
coeff = 2.0 * cosine
q0 = 0
q1 = 0
q2 = 0
for i in range(0, (int)(numSamples)):
#print("Hello")
q0 = (coeff * q1) - q2 + data[i]
q2 = q1
q1 = q0
#real = (q1 - (q2 * cosine)) / scalingFactor
#imag = (q2 * sine) / scalingFactor
#magnitude = math.sqrt((real * real) + (imag * imag))
magnitude = (q2*q2) + (q1*q1) - (coeff * q1 * q2)
return magnitude
max_mag_hi_freq = 0
maxvalue_mag_hi_freq = 0
for freq in hi_freqs:
mag = goertzel_mag(chunk_sample_count, freq, sample_rate, chunk_data)
#if c == 277:
# print("(%d)CHUNK %d magnitude with higher frequency (%d) = %f" % \
# (max_mag_hi_freq, c, freq, mag))
if(mag > maxvalue_mag_hi_freq):
maxvalue_mag_hi_freq = mag
max_mag_hi_freq = freq
max_mag_lo_freq = 0
maxvalue_mag_lo_freq = 0
for freq in lo_freqs:
mag = goertzel_mag(chunk_sample_count, freq, sample_rate, chunk_data)
#print("CHUNK %d magnitude with lower frequency (%d) = %f" % \
# (c, freq, mag))
if(mag > maxvalue_mag_lo_freq):
maxvalue_mag_lo_freq = mag
max_mag_lo_freq = freq
參考原文：https://blog.csdn.net/hankhanti/article/details/49902441
我們從1209, 1336, 1477和1633四個頻率對應的能量P中取最大值，記作Px，從679，770，852和941四個頻率對應的能量P中取出最大值Py。那麼Px和Py對應的頻率組合極有可能代表識別出一個DTMF符號。但是，我們還需要做一系列的判斷，來進一步評估：
- Px和Py是否足夠強大？我們可以設定一個門限，如果麼Px和Py這兩個任何一個低於這個門限，那麼N個取樣被評估為沒有識別出DTMF符號。參考資料[2]中建議這個門限值為4*105。但是如果取樣值的取值範圍是-32768到32767的話，實際上計算出來的P值會非常大，這個門限設為4*109都可以。
- Px和Py的差別是否太大？正常的DTMF訊號，這兩個能量應該接近，那麼如果差別較大，我們視為無效。參考資料[2]中建議的方法為：如果Py < Px * 0.398，那麼認為無效。如果Px < Py * 0.158也認為無效。但是實際上，我們將0.158改為0.5，識別效果更佳。
- 其它頻率的能量P有沒有很多接近Px和Py的？參考資料[2]中建議的方法為：首先取近Px和Py中較大的那個，設為Pm，如果其他頻率的P值有2個以上達到了Pm的15.8%，那麼認為是噪音導致，視為無效。
如果上述三個檢驗關卡都通過了，那麼我們可以將這N個取樣評估為包含一個DTMF符號，即Px和Py對應的頻率組合對應的某個符號。
參考資料：
解碼演算法：
/* Twist check
* CEPT => twist < 6dB
* AT&T => forward twist < 4dB and reverse twist < 8dB
* -ndB < 10 log10( v1 / v2 ), where v1 < v2
* -4dB < 10 log10( v1 / v2 )
* -0.4 < log10( v1 / v2 )
* 0.398 < v1 / v2
* 0.398 * v2 < v1
*/ if ( r[col] > r[row] )
{
/* Normal twist */
max_index = col;
if ( r[row] < (r[col] * 0.398) ) /* twist > 4dB, error */
see_digit = FALSE;
}
else /* if ( r[row] > r[col] ) */
{
/* Reverse twist */
max_index = row;
if ( r[col] < (r[row] * 0.158) ) /* twist > 8db, error */
see_digit = FALSE;
}