离散傅里叶变换与相关性计算

离散傅里叶变换（DFT）与相关性计算

本文主要介绍通过离散傅里叶变换的方法计算两个时间序列的相关性，不仅给出数学公式的推导，同时给出对应的 Matlab程序的验证。文中涉及到相关，卷积等基本概念，本文假设读者已经具备相关的数学基础。

更改说明：由于CSDN的MD语法在解析部分公式时出现错误，所以讲部分公式改为图片。
补充：关于循环卷积与频域的相乘关系

基本定义

首先，我们给出离散傅里叶变换，卷积与相关的定义：

离散傅里叶变换

长度为N的离散时间序列$x(n)$的离散傅里叶变换的定义为：
$$X(k)=\sum _{n=0}^{N-1}x(n)e^{-j2\pi nk/N}$$
对应的Matlab函数为fft(x)

卷积

长度为N的离散时间序列$x(n)$与$y(n)$的卷积定义为：
$$z(n)=\sum _{m=0}^{N-1}x(m)y(n-m)$$
得到的卷积$z(n)$的长度为$2N-1$，n的取值范围是0到2N-2，卷积对应的Matlab函数为conv(x, y)。举个简单的例子：

x = [1, 2, 3]
y = [1, -1]
z = conv(x, y) %  [1, 1, 1, -3]

相关

长度为N的离散时间序列$x(n)$的相关定义为：
$$z(n)=\sum _{m=0}^{N-1}x(m)y(m-n)$$
同样，得到的相关$z(n)$的长度为$2N-1$，卷积对应的Matlab函数为xcorr(x, y)， 注意n的范围是$-(N-1)$到$(N-1)$。举个简单的例子：

x = [1, 2, 3]
y = [1, -1, 1]
z = xcorr(x, y) %  [1, 1, 2, -1, 3]，对应的n是[-2, -1, 0, 1, 2]

注意xcorr中要求两个时间序列长度相等，如果不相等，会自动将较短的时间序列补0。

时域的卷积对应频域的相乘

设$z(n)$是时间序列$x(n)$与$y(n)$的卷积，长度为$2N-1$，则对应的离散傅里叶变换为：

$$\begin{array}{rl}Z(k)& =\sum _{n=0}^{2N-2}z(n)e^{-j2\pi nk/2N-1}\\ & =\sum _{n=0}^{2N-2}\sum _{m=0}^{N-1}x(m)y(n-m)e^{-j2\pi nk/2N-1}\end{array}$$
交换求和顺序得
$$\begin{array}{rl}Z(k)& =\sum _{m=0}^{N-1}x(m)\sum _{n=0}^{2N-2}y(n-m)e^{-j2\pi nk/2N-1}\\ & =\sum _{m=0}^{N-1}x(m)e^{-j2\pi mk/2N-1}\sum _{n=0}^{2N-2}y(n-m)e^{-j2\pi (n-m)k/2N-1}\\ & =\sum _{m=0}^{N-1}x(m)e^{-j2\pi mk/2N-1}\sum _{n=0}^{N-1}y(n)e^{-j2\pi nk/2N-1}\end{array}$$
注意等式（7）中第二项的变量替换。由于$Z(k)$序列的长度为$2N-1$，而$X(k)$与$Y(k)$的长度都为$N$，为了长度上一致，我们将$x(n)$与$y(n)$分别补0，使得最终的序列长度为$2N-1$，分别记作$x_p(n)$与$y_p(n)$，以$x_p(n)$为例对应的数学描述为:

$$x_p(n)=\{\begin{array}{rl}x(n)& n\le N-1\\ 0& n>N-1\end{array}$$
则其对应的傅里叶变换为：
$$X_p(k)=\sum _{n=0}^{N-1}{x}_p(n)e^{-j2\pi nk/2N-1}$$
联系等式(8)可以得到
$$Z(k)=X_p(k)Y_p(k)$$
我们使用一个简单的Matlab程序进行验证。

N = 20;
x = randn(1, N);
y = randn(1, N);

z = conv(x, y);

xp = [x, zeros(1, N-1)]; %补0
yp = [y, zeros(1, N-1)]; %补0

zk = fft(z);

xpk = fft(xp);
ypk = fft(yp);
zke = xpk.*ypk;%频域相乘

real_err = max(abs(real(zk)-real(zke)));
imag_err = max(abs(imag(zk) - imag(zke)));

fprintf('实补误差 %e\n虚部误差 %e\n',real_err, imag_err);

%% 程序输出结果
% 实部误差 1.421085e-14
% 虚部误差 8.881784e-15

其他的文献补0的个数都选择N而非N-1，背后的道理是一样的。选择N的好处是$X_p(k)$与$X(k)$有一个一致的对应关系，即$X_p(2k)=X(k)$

上面的定理可以简单描述为时域的卷积对应频域的乘积，一个典型应用加速卷积的计算。我们知道，如果直接套用卷积的定义，那么时间复杂度是$O(n^2)$，而快速傅里叶变换的时间复杂度是$O(n\log n)$，一个直接的思路就是利用傅里叶变换计算卷积，即首先根据$X_p(k)$与$Y_p(k)$得到时间$Z(k)$，然后进行傅里叶反变换得到卷积$z(n)$。

循环卷积与相乘的关系

在图像处理中，图像与卷积核的运算通常需要保持图像大小不变，这种情况下时域的卷积与频域的乘积的对应关系是否仍然保持呢？
Matlab对应的卷积函数为imfilter(u, v,'circular', 'conv')，给出简单的Matlab函数进行验证。

a=[1,2,3,4,5,6,7];
b = [1,2,3];
d = imfilter(a,b,'circular','conv');
% 结果 d => [ 11     8    12    16    20    24    21]

了解了循环卷积的运算过程，我们进行公式的推导。首先，对循环卷积进行表示。对于长度为N的序列$x(n)$其中，n=0, 1, 2, …N-1。以及长度为M的序列$h(n)$，这里对于$h(n)$，如果长度为奇数2M+1，则对应的n的取值范围为-M, -(M-1), … 0, …, (M-1), M。如果长度为偶数2M，则对应的n的范围表示为-M, …M-1。则循环卷积y(n)可以表示为：
$$y(n)=\sum _{m=0}^{N-1}x(m)h)(n-m),n=0,1,2,...N-1$$
根据上述公式，进行形同推导，则可以得到：
$$\begin{array}{rl}Y(k)& =\sum _{m=0}^{N-1}x(m)\sum _{n=0}^{N-1}h(n-m)e^{-j2\pi nk/N}\\ & =\sum _{m=0}^{N-1}x(m)e^{-j2\pi mk/N}\sum _{n=0}^{N-1}h(n-m)e^{-j2\pi (n-m)k/N}\end{array}$$
定义$h_p(n)$为将$h(n)$补0到长度为N，然后将向左循环平移M位。即：
$$h_p(n)=\{\begin{array}{rl}h(n)& 0\le n\le M-1\\ x(n-m+N)& n\ge n-M\\ 0& otherwise\end{array}$$
则$Y(k)=X(k)H_p(k)$，我们通过一个Matlab程序进行验证。

clc
clear
x=[1,2,3,4,5,6,7];
h = [1,2,-1];

y = imfilter(x, h, 'circular', 'conv'); %-3     6     8    10    12    14     9


hp = [h, zeros(1, length(x)-length(h))]; %补0
m = floor(length(h)/2);
hp = circshift(hp, -m); %向左循环平移
Yk = fft(x).*fft(hp);

ye = ifft(Yk);  %  -3.0000    6.0000    8.0000   10.0000   12.0000   14.0000    9.0000
% 计算结果可知，y与ye的结果一致(忽略舍入误差)

循环平移的傅里叶变换

循环平移的本质是将序列$x(n)$看做周期为N的周期序列，也就是N到2N-1的取值与0到N-1的取值有一一对应的关系。这里我们给出证明：

将序列$x(n)$向右做循环平移m位，得到y(n)，则：
$$y(n)=\{\begin{array}{rl}x(n-m)& n\ge m\\ x(n-m+N)& n<m\end{array}$$
则对应的傅里叶变换为：
$$\begin{array}{rl}Y(k)& =\sum _{n=0}^{N-1}y(n)e^{-j2\pi nk/N}\\ & =\sum _{n=m}^{N-1}x(n-m)e^{-j2\pi nk/N}+\sum _{n=0}^{m-1}x(n-m+N)e^{-j2\pi nk/N}\\ & =\sum _{n=m}^{N-1}x(n-m)e^{-j2\pi (n-m)k/N}{e}^{-j2\pi mk/N}+\sum _{n=0}^{m-1}x(n-m+N)e^{-j2\pi (n-m+N)k/N}{e}^{-j2\pi (m-N)k/N}\\ & =\sum _{n=0}^{N-1}x(n)e^{-j2\pi nk/N}{e}^{-j2\pi mk/N}\\ & =X(k)e^{-j2\pi mk/N}\end{array}$$

相关与卷积

熟悉卷积的读者应该知道卷积的计算过程分为两步：翻转和平移，与相关的计算相比，只是多了一步翻转操作。如果直接序列进行翻转，是否能用卷积代替相关操作呢？答案是是显然的，我们用一个Matlab程序做证明。

x = [1,2,1];
y = [1,2,3];

z1 = xcorr(x, y); %相关
z2 = conv(x, flip(y)); %先翻转，然后卷积

z1 % => [3, 8, 8, 4, 1]
z2 % => [3, 8, 8, 4, 1]

这里不做数学的证明，有兴趣的读者可以自己证明。证明的过程中注意相关计算结果的第一项对应的n是$-(N-1)$，而卷积操作的第一项对应的n是0。

通过傅里叶变换计算相关

如果直接利用相关的公式，时间复杂度是$O(n^2)$，而快速傅里叶变换的时间复杂度是$O(n\log n)$，那么能否使用fft方法来计算相关，从而提高运算速度呢？上一节我们已经知道相关与翻转后做卷积是等价的，那么一个直接的思路是：

1. 将序列$y(n)$翻转，得到$u(n)$，然后补0，得到$u_p(n)$，进行傅里叶变换得到$U_p(k)$
2. 将序列$x(n)$补0，得到$x_p(n)$，进行傅里叶变换得到 $X_p(k)$
3. 计算$Z(k)=X_{(}p)U_p(k)$，通过傅里叶反变换得到$z(n)$

上述计算的缺点是多了一个中间变量$u(n)$，那么能否用 $Y_p(k)$ 来计算$U_p(k)$呢？我们进行一个简单的推导。
$$\begin{array}{rl}{U}_p(k)& =\sum _{n=0}^{N-1}u(n)e^{-j2\pi nk/2N-1}\\ & =\sum _{n=0}^{N-1}y(N-1-n)e^{-j2\pi nk/2N-1}\\ & =\sum _{n=0}^{N-1}y(N-1-n)e^{-j2\pi (n+1-N)k/2N-1}{e}^{-j2\pi (N-1)k/2N-1}\\ & ={\left\{\sum _{n=0}^{N-1}y(N-1-n)e^{-j2\pi (N-1-n)k/2N-1}\right\}}^{\ast }{e}^{-j2\pi (N-1)k/2N-1}\\ & =Y_p^{\ast }(k)e^{-j2\pi (N-1)k/2N-1}\end{array}$$
其中$\ast$表示共轭运算，对应的matlab函数为conj(X)。

也就是说如果令$Z_1(k)=X_p(k)Y_p^{\ast }(k)$，对应的傅里叶反变换为$z_1(n)$，则$Z(k)=Z_1(k)e^{-j2\pi (N-1)k/2N-1}$，根据循环平移的傅里叶变换易知 $z(n)=z1(n-(N-1))$

因此，通过傅里叶变换计算相关运算的思路为：

1. 将$x(n)$与$y(n)$分别补N-1个0，得到$x_p(n)$与$y_p(n)$，进行快速傅里叶变换得到$X_p(k)$与$Y_p(k)$；
2. 计算$Z_1(k)=X(k)Y^{\ast }(k)$，并进行傅里叶反变换得到$z_1(n)$
3. 进行循环平移得到$z(n)$

同之前讨论过的一样，补0的个数只要大于等于N-1不影响计算结果，很多资料中都补N个0.

最后我们通过一个Matlab程序对比一下4种计算相关性的方法。

% 通过4种方式计算相关
%% 数据生成
N = 30;
x = randn(1, N);
y = [randn(1, 4), x(1:N-4)];
%y = randn(1, N);

%% 方法1， 通过Matlab自带函数xcorr进行计算

z1 = xcorr(x, y);

%% 方法2，通过翻转，求卷积

z2 = conv(x, flip(y));

%% 方法3，通过翻转，傅里叶变换计算
u = flip(y);
xp = [x, zeros(1, N-1)];
up = [u, zeros(1, N-1)];

zk = fft(xp).*fft(up);
z3 = ifft(zk);

%% 方法4，傅里叶变换方法
xp = [x, zeros(1, N-1)];

yp = [y, zeros(1, N-1)];

z1k = fft(xp).*conj(fft(yp));
z1n =  ifft(z1k);
z4 = circshift(z1n, N-1); %需要平移操作

%% 计算结果的可视化
lags = -(N-1):N-1;


scalediff = 5;
plot(lags, z1);
hold on;
plot(lags, z2 + scalediff*1);
plot(lags, z3 + scalediff*2);
plot(lags, z4 + scalediff*3);

plot([-4, -4], ylim, 'r:','LineWidth', 1.5); % 平移量
hold off
legend('方法一','方法二','方法三', '方法四');

结果：