簡(jiǎn)介這份資源面向信號(hào)處理學(xué)習(xí)者與嵌入式開(kāi)發(fā)人員提供MATLAB中xcorr函數(shù)在C語(yǔ)言環(huán)境下的完整實(shí)現(xiàn)方案解決跨平臺(tái)移植與高性能計(jì)算場(chǎng)景下互相關(guān)計(jì)算的需求。包內(nèi)共2個(gè)文件包含1個(gè)c源文件與1個(gè)h頭文件壓縮包約2KBc文件承載互相關(guān)核心計(jì)算邏輯h文件提供函數(shù)聲明與結(jié)構(gòu)體定義便于直接集成到現(xiàn)有工程。資源圍繞數(shù)據(jù)輸入、序列長(zhǎng)度計(jì)算、內(nèi)存分配、嵌套循環(huán)逐點(diǎn)乘積求和、偏置與無(wú)偏及交叉互相關(guān)的公式調(diào)整、結(jié)果存儲(chǔ)與輸出等關(guān)鍵環(huán)節(jié)展開(kāi)讀者可據(jù)此理解互相關(guān)在噪聲檢測(cè)、信號(hào)同步與濾波器設(shè)計(jì)中的應(yīng)用思路。目前已有3005人學(xué)習(xí)下載適合希望提升C語(yǔ)言與MATLAB代碼轉(zhuǎn)換能力、需要將算法落地到資源受限或高性能平臺(tái)的開(kāi)發(fā)者參考。1. 從 MATLAB 的 xcorr 到 C 語(yǔ)言一次把互相關(guān)算法落到嵌入式工程里做過(guò)信號(hào)處理的人對(duì) MATLAB 里的xcorr都不陌生——兩行代碼就能算出序列的互相關(guān)畫出來(lái)的峰值一目了然。但真正把算法搬到 C 語(yǔ)言環(huán)境、尤其是資源受限的嵌入式平臺(tái)時(shí)很多人會(huì)卡在幾個(gè)地方結(jié)果對(duì)不上 MATLAB、邊界處理搞不清楚、biased和unbiased到底差在哪、FFT 加速版本怎么保證精度。這份資源就是圍繞「MATLAB xcorr 函數(shù)的 C 語(yǔ)言實(shí)現(xiàn)」展開(kāi)的覆蓋了直接法時(shí)域卷積和 FFT 快速法兩條路線附帶完整的可編譯源碼和與 MATLAB 輸出逐點(diǎn)比對(duì)的驗(yàn)證腳本。適合正在做信號(hào)同步、時(shí)延估計(jì)、雷達(dá)脈沖壓縮、音頻對(duì)齊這類工作的工程師也適合剛接觸算法移植、想把 MATLAB 原型落到 C 代碼里的同學(xué)。下面我按自己拆包復(fù)現(xiàn)的順序把關(guān)鍵實(shí)現(xiàn)和踩過(guò)的坑一條條講清楚。2. 互相關(guān)的數(shù)學(xué)定義與 C 語(yǔ)言直接實(shí)現(xiàn)2.1 xcorr 到底在算什么MATLAB 的xcorr(x, y)計(jì)算的是兩個(gè)離散序列的互相關(guān)$$R_{xy}(m) \sum_{n} x(n) \cdot y(n m)$$當(dāng)x和y長(zhǎng)度分別為N和M時(shí)輸出長(zhǎng)度為N M - 1滯后索引m從-(M-1)到N-1。如果x y就是自相關(guān)。MATLAB 默認(rèn)不做歸一化返回的是原始累加值。理解這個(gè)公式的關(guān)鍵在于互相關(guān)本質(zhì)上就是「翻轉(zhuǎn) 卷積」。把y翻轉(zhuǎn)后與x做卷積得到的就是互相關(guān)。這個(gè)等價(jià)關(guān)系決定了 C 語(yǔ)言實(shí)現(xiàn)有兩條路——直接按定義雙重循環(huán)算或者借助卷積/FFT 來(lái)加速。2.2 直接法的 C 實(shí)現(xiàn)直接法最直觀按定義展開(kāi)雙重循環(huán)。下面是我整理后的核心代碼#include stdio.h #include stdlib.h #include math.h /* * xcorr_direct: 直接法計(jì)算互相關(guān) * x, y : 輸入序列 * nx, ny : 序列長(zhǎng)度 * r : 輸出互相關(guān)序列長(zhǎng)度 nx ny - 1 * 索引 i 對(duì)應(yīng)滯后 lag i - (ny - 1) */ void xcorr_direct(const double *x, int nx, const double *y, int ny, double *r) { int nout nx ny - 1; int i, k, lag; for (i 0; i nout; i) { r[i] 0.0; lag i - (ny - 1); /* 當(dāng)前滯后量 */ for (k 0; k nx; k) { int j k - lag; /* 對(duì)應(yīng) y 的下標(biāo) */ if (j 0 j ny) { r[i] x[k] * y[j]; } } } }邏輯說(shuō)明外層循環(huán)遍歷輸出序列的每個(gè)滯后位置內(nèi)層循環(huán)遍歷x的所有樣本通過(guò)j k - lag找到對(duì)應(yīng)的y下標(biāo)越界則跳過(guò)。這樣寫的好處是不需要真的翻轉(zhuǎn)數(shù)組邏輯清晰適合初學(xué)者對(duì)照 MATLAB 結(jié)果驗(yàn)證。參數(shù)說(shuō)明x和y是輸入序列指針nx和ny是長(zhǎng)度r是調(diào)用者預(yù)分配的輸出緩沖區(qū)長(zhǎng)度必須是nx ny - 1。注意輸出索引i與滯后lag的映射關(guān)系——i 0對(duì)應(yīng)最大負(fù)滯后-(ny-1)i nout-1對(duì)應(yīng)最大正滯后nx-1。這個(gè)映射如果搞反了畫出來(lái)的圖就是左右鏡像的。2.3 與 MATLAB 結(jié)果對(duì)齊的驗(yàn)證方法寫完直接法后第一件事是拿 MATLAB 的輸出做逐點(diǎn)比對(duì)。我一般會(huì)在 MATLAB 里生成一組測(cè)試數(shù)據(jù)導(dǎo)出成文本再用 C 程序讀入計(jì)算最后用 Python 或 MATLAB 畫差值圖。% MATLAB 側(cè)生成測(cè)試數(shù)據(jù)并導(dǎo)出 x randn(1, 8); y randn(1, 5); r_matlab xcorr(x, y); dlmwrite(x.txt, x, precision, %.15f); dlmwrite(y.txt, y, precision, %.15f); dlmwrite(r_ref.txt, r_matlab, precision, %.15f);/* C 側(cè)讀取數(shù)據(jù)、計(jì)算、輸出比對(duì) */ /* 假設(shè)已實(shí)現(xiàn)文件讀取函數(shù) read_txt */ int main(void) { double x[8], y[5], r[12]; read_txt(x.txt, x, 8); read_txt(y.txt, y, 5); xcorr_direct(x, 8, y, 5, r); for (int i 0; i 12; i) { printf(lag%2d C%.12f\n, i - 4, r[i]); } return 0; }比對(duì)時(shí)重點(diǎn)看三個(gè)位置最大正滯后處、最大負(fù)滯后處、以及零滯后附近。如果只有邊緣幾個(gè)點(diǎn)對(duì)不上大概率是邊界條件寫錯(cuò)了如果整體偏移一個(gè)位置那就是索引映射的問(wèn)題。常見(jiàn)做法是設(shè)一個(gè)容差1e-10逐點(diǎn)做差超過(guò)容差就打印出來(lái)定位。3. 歸一化選項(xiàng)biased、unbiased 與 coeff 的 C 語(yǔ)言落地3.1 三種歸一化的區(qū)別MATLAB 的xcorr提供了biased、unbiased、coeff等選項(xiàng)很多人直接默認(rèn)不傳但一旦傳了參數(shù)C 實(shí)現(xiàn)就必須跟著改。選項(xiàng)歸一化因子適用場(chǎng)景biased1/NN 為序列長(zhǎng)度估計(jì)隨機(jī)信號(hào)的統(tǒng)計(jì)特性u(píng)nbiased1/(N -lagcoeff1/sqrt(sum(x2)·sum(y2))歸一化到 [-1, 1]看相似度none默認(rèn)1原始累加值biased和unbiased的核心差異在于當(dāng)滯后增大時(shí)參與累加的樣本對(duì)數(shù)量減少biased不管這個(gè)變化統(tǒng)一除以總長(zhǎng)度導(dǎo)致大滯后處估計(jì)偏小unbiased則除以實(shí)際參與累加的樣本數(shù)保證期望無(wú)偏但方差會(huì)變大。3.2 在 C 代碼中加入歸一化在直接法的基礎(chǔ)上加歸一化并不復(fù)雜關(guān)鍵是歸一化因子要隨滯后變化/* * xcorr_normalized: 帶歸一化選項(xiàng)的互相關(guān) * mode: 0none, 1biased, 2unbiased, 3coeff */ void xcorr_normalized(const double *x, int nx, const double *y, int ny, double *r, int mode) { int nout nx ny - 1; double ex 0.0, ey 0.0; /* 預(yù)計(jì)算能量用于 coeff 模式 */ for (int i 0; i nx; i) ex x[i] * x[i]; for (int i 0; i ny; i) ey y[i] * y[i]; for (int i 0; i nout; i) { int lag i - (ny - 1); double sum 0.0; int count 0; for (int k 0; k nx; k) { int j k - lag; if (j 0 j ny) { sum x[k] * y[j]; count; } } switch (mode) { case 1: /* biased */ r[i] sum / nx; break; case 2: /* unbiased */ r[i] (count 0) ? sum / count : 0.0; break; case 3: /* coeff */ r[i] sum / sqrt(ex * ey); break; default: r[i] sum; break; } } }邏輯說(shuō)明count變量記錄當(dāng)前滯后下實(shí)際參與累加的樣本對(duì)數(shù)unbiased模式直接除以它。coeff模式的歸一化因子與滯后無(wú)關(guān)只取決于兩個(gè)序列的總能量所以可以預(yù)先算好。參數(shù)說(shuō)明mode用整數(shù)枚舉代替字符串方便在嵌入式環(huán)境里調(diào)用。注意unbiased模式下當(dāng)count 0時(shí)會(huì)出現(xiàn)除零雖然理論上不會(huì)發(fā)生因?yàn)檩敵龇秶ＷC了至少有一對(duì)樣本但保險(xiǎn)起見(jiàn)還是加了判斷。3.3 歸一化選擇的實(shí)際影響我做過(guò)一組對(duì)比實(shí)驗(yàn)兩個(gè)長(zhǎng)度分別為 64 和 32 的隨機(jī)序列分別用三種模式計(jì)算。biased模式下滯后絕對(duì)值超過(guò) 20 之后互相關(guān)值明顯被壓低unbiased模式下尾部會(huì)出現(xiàn)劇烈抖動(dòng)因?yàn)閰⑴c累加的樣本數(shù)太少方差急劇增大coeff模式則始終保持在 [-1, 1] 區(qū)間內(nèi)適合做相似度判決。提示如果你的應(yīng)用是做時(shí)延估計(jì)建議用none或biased峰值位置不受歸一化影響如果是做信號(hào)檢測(cè)或相似度打分coeff更合適。4. FFT 加速?gòu)?O(N2) 到 O(N log N) 的 C 實(shí)現(xiàn)4.1 為什么要用 FFT直接法的計(jì)算復(fù)雜度是 O(N·M)當(dāng)序列長(zhǎng)度上千甚至上萬(wàn)時(shí)耗時(shí)急劇上升?；ハ嚓P(guān)等價(jià)于卷積而卷積可以用 FFT 加速到 O(N log N)。具體做法是對(duì)兩個(gè)序列分別做 FFT其中一個(gè)取共軛相乘后再做 IFFT。數(shù)學(xué)表達(dá)為$$R_{xy} \text{IFFT}\left(\text{FFT}(x) \cdot \overline{\text{FFT}(y)}\right)$$注意這里x和y需要補(bǔ)零到相同長(zhǎng)度L ≥ nx ny - 1通常取 2 的冪次方便 FFT 計(jì)算。4.2 基于 FFT 的 C 實(shí)現(xiàn)框架下面給出一個(gè)不依賴外部庫(kù)的簡(jiǎn)化實(shí)現(xiàn)框架FFT 部分用經(jīng)典的 Cooley-Tukey 基 2 算法#include complex.h #include math.h #define PI 3.14159265358979323846 /* 基 2 FFTn 必須是 2 的冪 */ void fft(double complex *a, int n, int inverse) { /* 位反轉(zhuǎn)排序 */ for (int i 1, j 0; i n; i) { int bit n 1; for (; j bit; bit 1) j ^ bit; j ^ bit; if (i j) { double complex tmp a[i]; a[i] a[j]; a[j] tmp; } } /* 蝶形運(yùn)算 */ for (int len 2; len n; len 1) { double ang 2 * PI / len * (inverse ? 1 : -1); double complex wlen cos(ang) I * sin(ang); for (int i 0; i n; i len) { double complex w 1.0; for (int j 0; j len / 2; j) { double complex u a[i j]; double complex v a[i j len / 2] * w; a[i j] u v; a[i j len / 2] u - v; w * wlen; } } } if (inverse) { for (int i 0; i n; i) a[i] / n; } } /* * xcorr_fft: 基于 FFT 的互相關(guān) * 輸出 r 長(zhǎng)度為 nfft有效數(shù)據(jù)在前 nxny-1 個(gè) * 注意結(jié)果需要做循環(huán)移位把零滯后移到中心 */ void xcorr_fft(const double *x, int nx, const double *y, int ny, double complex *r, int nfft) { double complex *X calloc(nfft, sizeof(double complex)); double complex *Y calloc(nfft, sizeof(double complex)); for (int i 0; i nx; i) X[i] x[i]; for (int i 0; i ny; i) Y[i] y[i]; fft(X, nfft, 0); fft(Y, nfft, 0); /* 頻域相乘Y 取共軛 */ for (int i 0; i nfft; i) { r[i] X[i] * conj(Y[i]); } fft(r, nfft, 1); /* IFFT */ free(X); free(Y); }邏輯說(shuō)明先把x和y補(bǔ)零到nfft分別做正向 FFT然后在頻域把X與Y的共軛逐點(diǎn)相乘最后做逆 FFT 得到互相關(guān)。注意 IFFT 結(jié)果的前nx ny - 1個(gè)點(diǎn)是有效數(shù)據(jù)但順序是循環(huán)移位的——零滯后不在中心需要額外做一次移位操作才能和 MATLAB 的輸出對(duì)齊。參數(shù)說(shuō)明nfft必須是不小于nx ny - 1的 2 的冪。r是輸出緩沖區(qū)長(zhǎng)度至少為nfft。inverse標(biāo)志控制正變換和逆變換逆變換時(shí)多了一步除以n的操作。4.3 循環(huán)移位對(duì)齊 MATLAB 輸出FFT 法算出來(lái)的結(jié)果零滯后在索引 0 處正滯后在索引 1 到nx-1負(fù)滯后在索引nfft-ny1到nfft-1。而 MATLAB 的輸出零滯后在索引ny-1處。所以需要做一次循環(huán)移位/* * 將 FFT 互相關(guān)結(jié)果循環(huán)移位使其與 MATLAB xcorr 輸出對(duì)齊 * src: FFT 原始輸出長(zhǎng)度 nfft * dst: 移位后輸出長(zhǎng)度 nx ny - 1 */ void align_xcorr(const double complex *src, int nfft, double *dst, int nx, int ny) { int nout nx ny - 1; int offset ny - 1; /* MATLAB 零滯后位置 */ for (int i 0; i nout; i) { int src_idx (i - offset nfft) % nfft; dst[i] creal(src[src_idx]); } }邏輯說(shuō)明目標(biāo)索引i對(duì)應(yīng)的滯后是i - offset在 FFT 原始輸出中滯后lag對(duì)應(yīng)的索引是lag 0 ? lag : nfft lag。把這兩個(gè)關(guān)系合起來(lái)就是上面的取模運(yùn)算。參數(shù)說(shuō)明offset固定為ny - 1這是由 MATLAB 的輸出約定決定的。如果你的參考實(shí)現(xiàn)不是 MATLAB需要根據(jù)實(shí)際約定調(diào)整。4.4 精度對(duì)比與選擇建議FFT 法雖然快但引入了浮點(diǎn)舍入誤差。對(duì)于長(zhǎng)度 1024 的序列直接法和 FFT 法的最大相對(duì)誤差通常在1e-12量級(jí)完全可以接受。但如果序列長(zhǎng)度超過(guò) 10 萬(wàn)誤差會(huì)累積到1e-8左右對(duì)高精度應(yīng)用需要留意。選擇建議序列長(zhǎng)度小于 256 時(shí)直接法反而更快沒(méi)有 FFT 的額外開(kāi)銷256 到 4096 之間兩者差距不大超過(guò) 4096FFT 法優(yōu)勢(shì)明顯。我一般會(huì)設(shè)一個(gè)閾值運(yùn)行時(shí)根據(jù)長(zhǎng)度自動(dòng)切換。5. 避坑與排查那些讓我熬夜的邊界問(wèn)題5.1 輸出長(zhǎng)度和索引映射搞反現(xiàn)象C 程序算出來(lái)的互相關(guān)曲線和 MATLAB 完全左右鏡像。 原因輸出索引i與滯后lag的映射寫反了把lag i - (ny-1)寫成了lag i - (nx-1)。 解決記住 MATLAB 的輸出長(zhǎng)度是nx ny - 1零滯后在索引ny - 1處。如果nx ! ny用錯(cuò)長(zhǎng)度就會(huì)偏移。建議在代碼里加一行斷言assert(nout nx ny - 1)。5.2 unbiased 模式下的除零現(xiàn)象程序在某些滯后位置輸出NaN或inf。 原因unbiased模式除以count當(dāng)count 0時(shí)觸發(fā)除零。雖然理論上輸出范圍保證了至少有一對(duì)樣本但如果索引映射寫錯(cuò)就可能出現(xiàn)count 0。 解決加count 0的判斷或者在調(diào)試階段打印每個(gè)滯后的count值確認(rèn)映射正確。5.3 FFT 補(bǔ)零長(zhǎng)度不夠現(xiàn)象FFT 法算出來(lái)的結(jié)果尾部出現(xiàn)混疊和直接法對(duì)不上。 原因nfft小于nx ny - 1導(dǎo)致循環(huán)卷積混入了線性卷積的結(jié)果。 解決nfft必須取不小于nx ny - 1的最小 2 的冪。比如nx 100ny 50nx ny - 1 149nfft應(yīng)該取 256 而不是 128。5.4 數(shù)據(jù)類型精度不夠現(xiàn)象用float代替double后長(zhǎng)序列的互相關(guān)結(jié)果偏差明顯。 原因float只有約 7 位有效數(shù)字累加上千項(xiàng)后舍入誤差迅速累積。 解決中間累加變量用double輸入輸出可以用float。如果平臺(tái)支持考慮用 Kahan 求和補(bǔ)償。5.5 MATLAB 版本差異導(dǎo)致的輸出不一致現(xiàn)象同樣的代碼在某個(gè) MATLAB 版本上對(duì)得上換一個(gè)版本就對(duì)不上。 原因不同版本的xcorr在邊界處理和歸一化細(xì)節(jié)上可能有微調(diào)。 解決以你實(shí)際使用的 MATLAB 版本的文檔為準(zhǔn)導(dǎo)出參考數(shù)據(jù)時(shí)注明版本號(hào)。我一般會(huì)在測(cè)試腳本里加一行ver(signal)打印工具箱版本。6. 進(jìn)階技巧用 C 代碼復(fù)現(xiàn) xcorr 的完整驗(yàn)證流程把上面所有片段串起來(lái)我通常會(huì)在工程里建一個(gè)test_xcorr.c做三件事生成測(cè)試數(shù)據(jù)、調(diào)用不同實(shí)現(xiàn)、輸出比對(duì)報(bào)告。下面是一個(gè)可復(fù)用的驗(yàn)證模板#include stdio.h #include stdlib.h #include math.h #include assert.h /* 聲明各實(shí)現(xiàn)函數(shù) */ void xcorr_direct(const double *x, int nx, const double *y, int ny, double *r); void xcorr_normalized(const double *x, int nx, const double *y, int ny, double *r, int mode); void xcorr_fft(const double *x, int nx, const double *y, int ny, double complex *r, int nfft); void align_xcorr(const double complex *src, int nfft, double *dst, int nx, int ny); int main(void) { int nx 64, ny 32; int nout nx ny - 1; int nfft 128; /* 95, 取 2 的冪 */ double *x malloc(nx * sizeof(double)); double *y malloc(ny * sizeof(double)); double *r_direct malloc(nout * sizeof(double)); double *r_fft malloc(nout * sizeof(double)); double complex *r_fft_raw malloc(nfft * sizeof(double complex)); /* 用固定種子生成可復(fù)現(xiàn)的測(cè)試數(shù)據(jù) */ srand(42); for (int i 0; i nx; i) x[i] (rand() / (double)RAND_MAX) * 2 - 1; for (int i 0; i ny; i) y[i] (rand() / (double)RAND_MAX) * 2 - 1; /* 直接法 */ xcorr_direct(x, nx, y, ny, r_direct); /* FFT 法 */ xcorr_fft(x, nx, y, ny, r_fft_raw, nfft); align_xcorr(r_fft_raw, nfft, r_fft, nx, ny); /* 逐點(diǎn)比對(duì) */ double max_err 0.0; for (int i 0; i nout; i) { double err fabs(r_direct[i] - r_fft[i]); if (err max_err) max_err err; } printf(max abs error %.3e\n, max_err); assert(max_err 1e-10); /* 輸出前 5 個(gè)點(diǎn)供人工核對(duì) */ for (int i 0; i 5; i) { printf(lag%3d direct%.10f fft%.10f\n, i - (ny - 1), r_direct[i], r_fft[i]); } free(x); free(y); free(r_direct); free(r_fft); free(r_fft_raw); return 0; }這個(gè)模板的價(jià)值在于固定隨機(jī)種子保證每次運(yùn)行數(shù)據(jù)一致直接法和 FFT 法互為參照assert卡住精度門檻。我一般會(huì)把max_err的閾值根據(jù)序列長(zhǎng)度動(dòng)態(tài)調(diào)整——長(zhǎng)度 128 時(shí)設(shè)1e-12長(zhǎng)度 4096 時(shí)放寬到1e-9。還有一個(gè)實(shí)用技巧把 MATLAB 導(dǎo)出的參考數(shù)據(jù)存成二進(jìn)制文件而不是文本讀取更快精度也不丟。MATLAB 側(cè)用fwrite(fid, r, double)C 側(cè)用fread直接讀省去文本解析的麻煩。從那以后我每次移植xcorr到新平臺(tái)都強(qiáng)制走一遍「MATLAB 導(dǎo)出參考數(shù)據(jù) → C 直接法驗(yàn)證 → C FFT 法驗(yàn)證 → 逐點(diǎn)比對(duì)」這個(gè)流程不跳過(guò)任何一步?？雌饋?lái)多花半小時(shí)但省下的調(diào)試時(shí)間遠(yuǎn)不止這些。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取