MATLAB中circshift函数的原理分析——psf2otf函数的核心

之所以讲到MATLAB中circshift函数，也是源于Rafael Gonzalez的这个图，作为前几篇答廖老师问的blog的基础。
Rafael Gonzalez的这个图无论从哪幅图到哪幅图都不是直接的傅里叶变换或傅里叶逆变换，需要循环移位，即circshift函数。

这就需要从头说起。

离散时间傅里叶变换（DTFT）的一个性质是，如果一个序列共轭对称，即$x[n]=x^{\ast }[-n]$，那么它的傅里叶变换是实数。因此，当fft函数的输出结果出乎意料地为复数时，人们有时会感到惊讶。

用实序列说明。当输入序列是实偶对称时，其傅里叶变换的输出序列也是实数且偶对称的。

例如：

n = -3:3;
x = exp(-abs(n)/2)x =0.2231    0.3679    0.6065    1.0000    0.6065    0.3679    0.2231

这个序列看起来是对称的，但当我们计算fft(x)时，其输出却是复数。

>> X = fft(x)X =列 1 至 53.3951 + 0.0000i  -1.0726 - 0.5166i   0.2154 + 0.2701i  -0.0593 - 0.2598i  -0.0593 + 0.2598i列 6 至 70.2154 - 0.2701i  -1.0726 + 0.5166i

原因是fft函数计算的是在区间$0\le n<N$内非零的序列$x[n]$的离散傅里叶变换。实际上，$x[n]$并不是关于原点对称的，它是一个对称序列的移位版本，而这种移位导致了fft的输出为复数。

为了得到实数值的傅里叶变换，我们需要将序列循环移位，使中心元素移动到向量的左边。对于长度为奇数或偶数的序列，我们都可以使用以下方法来找到中心位置，并进行相应的循环移位：

xs = circshift(x, [0 -floor(length(x)/2)]);
xs =1.0000    0.6065    0.3679    0.2231    0.2231    0.3679    0.6065

现在，如果我们计算 fft(xs)，将会得到预期的实数输出。

>> Xs = fft(xs)Xs =3.3951    1.1905    0.3454    0.2665    0.2665    0.3454    1.1905

结合之前讲的循环移位，用这幅图解释，循环移位后是最后的结果，但是对应于第三幅中中心在原点的情况，这才是偶函数。

如果打算进行零填充（zero-padding），应该先填充再应用循环移位（这里很关键，我曾经在这里错过，百思不得其解，思考了数年之久）。下面是一个示例，其中我们将原始序列扩展到长度为128，并对其进行处理：

x128 = x;
x128(128) = 0; % 零填充到长度128
x128s = circshift(x128, [0 -floor(length(x)/2)]); % 循环移位
X128 = fft(x128s);

检查是否为实数：

isreal(X128)ans =logical0

这里可能会出现一个小的误差，认为这个过程不能产生实数结果，但实际上是因为浮点舍入误差。查看虚部的大小可以发现它们非常小：

max(abs(imag(X128(:))))ans =2.1252e-16

可以使用real函数去除可以忽略不计的虚部。下面是绘制实值傅里叶变换的方法。我将使用频率轴标记技术来展示结果：

w = unwrap(fftshift(2*pi * (0:(128-1)) / 128) - 2*pi);
figure; plot(w/pi, fftshift(real(X128)), "LineWidth", 1);
xlabel('弧度 / \pi');
box off;
grid on;