Periodic Signal

演算法（Fourier Analysis）

推廣到複數是Lapalce Analysis。

演算法（Wavelet Analysis）

Fourier transform：波形分解成一群簡諧波，頻率為整數倍。

傅立葉轉換的計算結果，經常轉換成頻譜。除了廣為人知的頻譜之外，還有許多衍生產物，命名方式是翻轉字母順序，非常搞笑。

specmurt和cepstrum主要用於共鳴聲音，例如樂器聲、說話聲。共鳴產生諧音：強度頻譜上，強度規律地出現，其頻率呈倍數。頻率軸取log，讓出現間隔成為等距，叫做specmurt。強度規律地出現，宛如波，於是有人再度套用一次傅立葉轉換找到基頻（個人認為莫名其妙），叫做cepstrum。由於頻譜有一些重大缺點，又加上計算時間長，所以並不好用。

spectrum：頻譜。分為強度頻譜、相位頻譜。
specmurt：名稱不詳。頻譜的座標軸，頻率取log，符合人類聽覺感受。
cepstrum：倒頻譜。強度頻譜，強度值取log，實施（逆向）傅立葉轉換。

const float π = 3.1415026; void ft(float xr[], float xi[], int n, float yr[], float yi[] ) { float ω = 2.0 * π / n; for (int i=0; i<n; i++) { yr[i] = 0; yi[i] = 0; for (int j=0; j<n; j++) { float er = cos(-ω*i*j); float ei = sin(-ω*i*j); yr[i] += xr[j] * er - xi[j] * ei; yi[i] += xr[j] * ei + xi[j] * er; } } } void ift(float xr[], float xi[], int n, float yr[], float yi[] ) { float ω = 2.0 * π / n; for (int i=0; i<n; i++) { yr[i] = 0; yi[i] = 0; for (int j=0; j<; j++) { float er = cos(ω*i*j); float ei = sin(ω*i*j); yr[i] += xr[j] * er - xi[j] * ei; yi[i] += xr[j] * ei + xi[j] * er; } } } void spectrum(float real[], float imag[], int n, float magnitude[], float phase[] ) { for (int i=0; i<n; ++i) { magnitude[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]) / n; phase[i] = atan2(imag[i], real[i]); } } void cepstrum(float magnitude[], int n, float cepstrum[] ) { float ar[n], ai[n], br[n], bi[n]; for (int i=0; i<n; ++i) ar[i] = log(magnitude[i]); for (int i=0; i<n; ++i) ai[i] = 0; ift(ar, ai, n, br, bi); for (int i=0; i<n; ++i) cepstrum[i] = sqrt(br[i]*br[i] + bi[i]*bi[i]); }

Peak Detection

找到波形的尖峰。

一、消除鋸齒：方法很多，諸如
　口、時域的平滑效果（k點平均數）（k點中位數）。
　口、頻域的刪除高頻（高頻形成鋸齒）。
　口、時域的Linear Prediction（迴歸函數）。
二、尋找極值：中央高、兩側低。
三、估計峰尖：極值前後，三點做拋物線內插，求出拋物線頂點。

void peak_detection(float src[], int n) { /* 5-point moving average */ const int k = 5; const int ½k = k / 2; float y[n]; float sum = 0; for (int i = -½k, p = -k, q = 0; i<n; ++i) { if (p >= 0) sum -= src[p]; if (q < n) sum += src[q]; p++; q++; if (i < 0) continue; if (i < ½k) y[i] = sum / q; else if (i >= n-½k) y[i] = sum / (n-p); else y[i] = sum / k; } /* find index of maximum */ vector<float> peak; for (int i=1; i<n-1; ++i) if (y[i] > y[i-1] && y[i] > y[i+1]) peak.push_back(i); /* parabola interpolation */ for (int i=0; i<peak.size(); ++i) { int j = peak[i]; if (src[j] > src[j-1] && src[j] > src[j+1]) peak[i] = (src[j+1] - src[j-1]) / 2.0f / (src[j] + src[j] - src[j+1] - src[j-1]); else // 當初smooth做的不夠，導致極值是凹的。 peak[i] = max(src[j-1], src[j+1]); } }

Peak Continuation

尖峰隨著時間連續變化。追蹤每個時刻的尖峰位置。

Frequency Detection

找到波形的頻率。

時域的方法
peak detection：找到兩個波峰，位置相減得到波長，波長倒數得到頻率。僅適合純音。
zero-crossing rate：波形穿越零的次數。僅適合純音。
 ACF：位移、相乘、加總（內積）。各種位移量，找最大值。
AMDF：位移、相減再絕對值、加總（絕對值誤差）。各種位移量，找最小值。
 YIN：位移、相減再平方、加總（平方誤差）。各種位移量，除以累積和，找最小值。
複合：例如 ACF / (AMDF + 1.0)，找最大值。

頻域的方法
peak detection：找到強度最高的頻率。計算速度較慢。

// zero-crossing rate float frequency_detection(float src[], int n) { int zero = 0; // for (int i=0; i<fs-1; ++i) for (int i=0; i<n-1; ++i) { if ((src[i] > 0 && src[i+1] <= 0) || (src[i] < 0 && src[i+1] >= 0)) zero++; } // return (float)zero / 2.0f; return (float)zero / n * fs / 2.0f; } // ACF & AMDF & YIN float frequency_detection(float src[], int n) { int gap = 0; float max = -1e18; float min = +1e18; float yin_sum = 0; for (int i=1; i<n/2; ++i) // gap { float acf = 0, amdf = 0, yin = 0; for (int j=0; j+i<n; ++j) { acf += pcm[j] * pcm[j+i]; amdf += fabs(pcm[j] - pcm[j+i]); yin += (pcm[j] - pcm[j+i]) * (pcm[j] - pcm[j+i]); } // if (acf > max) max = acf , gap = i; // if (amdf < min) min = amdf, gap = i; // yin_sum += yin; // float r = yin / yin_sum; // if (r < min) min = r , gap = i; float r = acf / (amdf + 1.0); if (r > max) max = r , gap = i; } return (float)fs / gap; }

Envelope Detection

找出波形的包絡線。

波形劇烈變動，處處尖峰；波形簡化為包絡線，容易辨認尖峰。

聲音訊號波形的包絡線（時域）：
　回、正向傅立葉轉換，數列後半（共軛對稱部分）設為0，
　　　逆向傅立葉轉換，取絕對值，乘以2。
　回、Hilbert transform，取絕對值。
　　　效果同上。

                                   |  Assume g(t) x cos(wt)
1. fourier transform               |  G(f) ⋅ 0.5 ⋅ (delta(w) + delta(-w))
2. negative frequency: value -> 0  |  G(f) ⋅ 0.5 ⋅ delta(w)
3. inverse fourier transform       |  0.5 g(t) x e𝑖wt
4. abs                             |  0.5 |g(t)|
5. multiply 2                      |  |g(t)|

void envelope(float frame[], int n, float env[]) { float ar[n], ai[n], br[n], bi[n]; for (int i=0; i<n; ++i) ai[i] = 0; fft(frame, ai, n, br, bi); for (int i=n/2; i<n; ++i) br[i] = bi[i] = 0; ifft(br, bi, n, ar, ai); for (int i=0; i<n; ++i) env[i] = sqrt(ar[i]*ar[i] + ai[i]*ai[i]) / n * 2; }

強度頻譜的包絡線（頻域）：
　回、聲音訊號實施linear prediction，轉換到頻域。稱作LPC spectrum。
　　　linear prediction的項數，設定成尖峰數量的兩倍（共軛對稱）。pole即尖峰。
　回、強度頻譜套用lowpass filter。因其具有平滑效果。
　　　不合邏輯，比較少用。
　回、承上，改成取log，並且改成在頻域實施lowpass filter。
　　　此即cepstrum乘上0或1。
　　　不合邏輯，比較少用。

https://www.phon.ucl.ac.uk/courses/spsci/dsp/lpc.html

Hibert–Huang transform：找出波形的包絡線的中間線。

Energy Detection

找出波形的能量。

Spectrum Sensing: Enhanced Energy Detection Technique Based on Noise Measurement

Noise（Random Signal）

隨機訊號，通常稱作「雜訊」或「噪訊」。

由於沒人研究隨機數列，大家只好援引浮動數列，以描述雜訊。

例如「White Gauss Noise」：高斯隨機雜訊，每個隨機變數的平均數（浮動中心）和變異數（浮動範圍）均相同。因為剛好是白雜訊，所以名稱裡面有白。

Noise的頻譜

數學家仿照光譜由紅到紫的特性，嘗試分類雜訊。

 white: 強度為常數
  grey: 強度符合人類聽覺曲線。（不那麼白）
   red: 強度正比於頻率倒數平方。　　以頻率對數為座標軸，漸減6dB。
  pink: 強度正比於頻率倒數。　　　　以頻率對數為座標軸，漸減3dB。（不那麼紅）
violet: 強度負正比於頻率倒數平方。　以頻率對數為座標軸，漸增6dB。
  blue: 強度負正比於頻率倒數。　　　以頻率對數為座標軸，漸增3dB。（不那麼紫）

http://vellocet.com/dsp/noise/VRand.html

Smooth Noise

計算學家運用「內插」，製造柔順的雜訊。網路上已有詳細教學文章，請讀者自行參考。

value noise：等距設置隨機數值，內插得到其餘數值。內插是為了製造綿延感。

perlin noise：等距設置遞迴碎形，內插得到其餘數值。有些規律、又不失隨機；好像有道理、又好像在鬼扯。非常奇葩。

http://mrl.nyu.edu/~perlin/noise/

fBm noise：疊加各種解析度（頻率）的雜訊。解析度是2的各種次方。雜訊可以是上述任意一種。隱含混沌與碎形的概念，更細膩、更自然。

wavelet noise：疊加各種頻率暨振幅的波。類似fBm noise。