第19章：用FFT计算线性卷积和循环互相关

陈老湿

中国科学技术大学信息与通信工程硕士

高西全《数字信号处理——原理、实现与应用》书籍中讲到 “ 频谱分析和滤波是最基本的信号处理，频谱分析就是计算信号的离散傅里叶变换（DFT），滤波实质上就是计算两个信号的卷积。 ”

FFT无疑是数字信号处理课程中的一个重点，最近碰到了用FFT计算两个序列的相关性，于是将FFT计算相关和卷积知识复习了下。

本章主要内容如下：

一、用FFT计算卷积

二、用FFT计算相关

三、推荐阅读与总结

一、用FFT计算线性卷积

假设h(n)的长度为N，x(n)的长度为M，则h(n)与x(n)线性卷积出来结果序列的长度为N+M-1。

线性卷积： \[y\left( n \right) = h\left( n \right) * x\left( n \right) = \sum\limits_{m = 0}^{N - 1} {h\left( m \right)x\left( {n - m} \right)} \]

循环卷积： \[{y_c}\left( n \right) = h\left( n \right) \otimes x\left( n \right) = \sum\limits_{m = 0}^{N - 1} {h\left( m \right)x{{\left( {\left( {n - m} \right)} \right)}_L}} {R_L}\left( n \right)\]

其中 \[{x{{\left( {\left( {n - m} \right)} \right)}_L}}\] 是 \[{x\left( {n - m} \right)}\] 以L为周期的周期延拓序列，注意有个矩形窗 \[{R_L}\left( n \right)\] 。 \[L \ge N + M - 1\] 是循环卷积结果和线性卷积结果相等的充要条件。

值得注意的是：循环卷积，一定要指明做L点的循环卷积；且 \[L \ge N + M - 1\] 这个条件求职笔试时，喜欢考填空题。

我目前仿真过程中，碰到用线性卷积的两个知识点：

（1）卷积码，这个可以参考《陈老湿：第2章（2）：卷积码》；

（2）信号经过滤波器，比如脉冲成型滤波器、低通滤波器，也是卷积过程，MATLAB中可以使用conv函数，或者filter函数，可以参考《陈老湿：第1章（1）：BPSK调制解调器仿真》。

%%%%%%%%%%%%%%%%%%%%%  计算线性卷积 %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  test_linear_conv.m  %%%%%%%%%
%%%%% date:2022年3月5日  author:飞蓬大将军 %%%%%%%%%%
%%%%%程序说明
%%%对比用conv和fft来计算卷积
A = [5 6 7 8 9 10]; 
B = [1 2 3 4] ;
length_total = length(A) +  length(B) +1; %L的长度，需要>=N+M-1
A_temp = [A zeros(1,length_total-length(A))];
B_temp = [B zeros(1,length_total-length(B))];
C_linear = conv(A,B); %使用conv函数
C_fft_temp = ifft(fft(A_temp).*(fft(B_temp))); %使用循环卷积计算线性卷积
C_fft = C_fft_temp(1,1:length(A) +  length(B) -1); %取主值
%%%结论
%**当L>=N+M-1，C_linear和C_fft两者结果才一致

二、用FFT计算循环互相关

用FFT计算相关，源于我之前看别人计算两个序列的互相关性，会使用下面这一种计算形式：

auto_corr_temp1 = ifft(fft(seq_polar_temp).*conj(fft(seq_polar)));

和循环卷积来计算线性卷积的代码稍有不同， 多了一个取共轭的操作，即conj 。

下面是MATLAB中对xcorr函数的介绍，不过下面这张截图中不是计算循环互相关。

A = [5 6 7 8]; 
B = [1 2 3 4] ;
C_fft = ifft(fft(A).*conj(fft(B))); %用fft实现
C_corr = zeros(1,length(A));
for iii = 1:length(A)
    A_temp = circshift(A,iii-1);
    C_corr(1,iii) = sum(A_temp.*B);
end

从结果上来看， C_corr = C_fft = [70 64 62 64]。详细过程手动如下：

C_corr(1,1) = 5*1+6*2+7*3+8*4 = 70；C_corr(1,2) = 5*2+6*3+7*4+8*1 = 64；

C_corr(1,3) = 5*3+6*4+7*1+8*2 = 62；C_corr(1,2) = 5*4+6*1+7*2+8*3 = 64。

那么计算循环互相关，有什么作用呢？

比如在导航信号当中，假设每一个比特都乘以相同的扩频码，如果接收机需要估计接收信号与本地产生的扩频码序列相差了多少码相位？

接收机便可以利用fft的方法计算循环互相关，并观察最大值的位置得知相差的码相位。 在导航信号中，之所以使用循环互相关，是因为每一个比特乘以相同的扩频码，且信号到达接收机的时刻是不确定的。

%%%%%%%%%%%%%%%%%%%%%  测试m序列的自相关性 %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%  test_m_sequence.m  %%%%%%%%%
%%%%% date:2022年3月5日  author:飞蓬大将军 %%%%%%%%%%
%%%%%程序说明
%%%用fft和cicrshift两种计算方式测试m序列相关性
%***************************************************
% 前半部分来自：https://blog.csdn.net/cjbct/article/details/78153616
% 此函数用来生成m序列
% coef为反馈系数向量
% Author: FastestSnail
% Date: 2017-10-03
%***************************************************
coef = [0 0 1 1];
m=length(coef);
len=2^m-1; % 得到最终生成的m序列的长度     
backQ=0; % 对应寄存器运算后的值，放在第一个寄存器
seq=zeros(1,len); % 给生成的m序列预分配
registers = [1 zeros(1, m-2) 1]; % 给寄存器分配初始结果
for i=1:len
    seq(i)=registers(m);
    backQ = mod(sum(coef.*registers) , 2); %特定寄存器的值进行异或运算，即相加后模2
    registers(2:length(registers)) = registers(1:length(registers)-1); % 移位
    registers(1)=backQ; % 把异或的值放在第一个寄存器的位置
%% 测试m序列循环互相关性
%*** 一般以本地信号为基准
seq_polar_local = 2*seq-1;  %假设是本地
location = 3;
seq_polar_signal = circshift(seq_polar_local,location); %假设外部信号进来
auto_corr_temp1 = ifft(fft(seq_polar_local).*conj(fft(seq_polar_signal)));
% auto_corr_temp1 = ifft(fft(seq_polar_signal).*conj(fft(seq_polar_local)));
index1 = find(auto_corr_temp1 == max(auto_corr_temp1))-1; %标号从1开始
auto_corr_temp2 = zeros(1,length(seq_polar_signal));
for iii = 1:length(seq_polar_signal)
    seq_temp = circshift(seq_polar_signal,iii-1);
    auto_corr_temp2(1,iii) = sum(seq_polar_local.*seq_temp);
index2 = find(auto_corr_temp2 == max(auto_corr_temp2))-1;%标号从1开始
stem(auto_corr_temp1);
%*** 以外部信号为基准
% for iii = 1:length(seq_polar_signal)
%     seq_temp = circshift(seq_polar_local,iii-1);
%     auto_corr_temp2(1,iii) = sum(seq_polar_signal.*seq_temp);
% end