将快速傅里叶变换(FFT)转换为定点

打开实时脚本

这个例子展示了如何将教科书版本的快速傅里叶变换(FFT)算法转换成定点MATLAB®代码。

运行以下代码以捕获当前状态，并重置全局状态。

FIPREF_STATE = get(fipref);重置(fipref)

教科书FFT算法

FFT是一个从时域到频域的复值线性变换。例如，如果你构造一个矢量作为两个正弦信号的和，并用FFT变换它，你可以在FFT幅值图中看到频率的峰值。

N = 64;%点数Fs = 4;%采样频率，单位为Hzt = (0:(n-1))/Fs;%时间矢量f = linspace(0,Fs,n);频率矢量F0 = 0.2;F1 = 0.5;%频率，单位为HzX0 = cos(2* *f0*t) + 0.55*cos(2* *f1*t);%时域信号X0 =复数(X0);教科书算法要求输入是复杂的Y0 = fft(x0);频域变换fft()是MATLAB的一个内置函数fi_fft_demo_ini_plot (t, x0, f, y0);绘制fft和时域信号的结果

图中包含2个轴对象。axis对象1包含一个类型为line的对象。该对象表示x0。Axes对象2包含一个类型为line的对象。该对象表示abs(fft(x0))。

频率图中0.2和0.5 Hz处的峰值对应于时域信号在该频率处的两个正弦曲线。

注意3.5和3.8 Hz的反射峰。当FFT的输入是实值时，就像在这个例子中一样，那么输出 $y$ 是共轭对称的:

$y （ k ）＝连词（ y （ n - k ））$

FFT有许多不同的实现，每一种都有自己的成本和收益。您可能会发现不同的算法比这里给出的算法更适合您的应用程序。该算法为您提供了一个如何开始自己探索的示例。

这个例子使用了在本书第45页的算法1.6.2中显示的实时抽取单位步幅FFT快速傅里叶变换的计算框架作者:查尔斯·范洛安

在伪代码中，教材中的算法如下:

1.6.2算法。如果 $x$ 是长度的复向量吗 $n$ 而且 $n ＝ 2^{t}$ ，则以下算法覆盖 $x$ 与 $F_{n} x$ ．

$x ＝ P_{n} x$

$w ＝ w_{n}^{（长）}$ (见Van Loan章节1.4.11)

为 $问＝ 1 ： t$

$l ＝ 2^{问} ； r ＝ n / l ； l_{＊} ＝ l / 2 ；$

为 $k ＝ 0 ： r - 1$

为 $j ＝ 0 ： l_{＊} - 1$

$τ ＝ w （ l_{＊} - 1 + j ） \cdot x （吉隆坡 + j + l_{＊} ）$

$x （吉隆坡 + j + l_{＊} ）＝ x （吉隆坡 + j ） - τ$

$x （吉隆坡 + j ）＝ x （吉隆坡 + j ） + τ$

结束

教科书算法使用从零开始的索引。 $F_{n}$ 是一个n × n的傅里叶变换矩阵， $P_{n}$ 是n × n位反转排列矩阵，和 $w$ 是旋转因子的复向量。旋转因素， $w$ ，为按以下算法计算的单位的复根:

函数W = fi_radix2twiddles(n)@ FI_RADIX2TWIDDLES基数-2 FFT示例的旋转因子。% W = FI_RADIX2TWIDDLES(N)计算长度为N-1的向量W%旋转因子在FI_M_RADIX2FFT示例代码中使用。％%参见FI_RADIX2FFT_DEMO。%的参考:％算法1.6.2的%旋转因子，第45页，Charles Van Loan，快速傅里叶变换的计算框架，SIAM，费城，1992年。％版权所有The MathWorks, Inc.％T = log2(n);如果地板(t) ~= t误差(N一定是2的幂）;结束W = 0 (n-1,1);K = 1;L = 2;公式1.4.11，第34页而L <= n = 2*pi/L;%算法1.4.1，第23页为j = 0: (L / 2 - 1) w (k) =复杂(cosθ(j *), sinθ(j *));K = K + 1;结束L = L*2;结束

Figure (gcf) CLF w0 = fi_radix2twiddles(n);polarplot(角(w0), abs (w0),“o”)标题(“旋转因素:统一的复根”）

图中包含一个axes对象。axis对象包含一个类型为line的对象。

验证浮点代码

要在MATLAB®中实现算法，您可以使用fi_bitreverse函数对输入序列进行位反转。必须向索引中添加1，才能将它们从从零开始转换为从一开始。

函数X = fi_m_radix2fft_algorithm1_6_2(X, w)%FI_M_RADIX2FFT_ALGORITHM1_6_2基数-2 FFT示例。的基数-2 fft_algorithm1_6_2 (X, W)计算的FFT%的输入向量X加上旋转因子w，假设输入X为%复杂。％向量X的长度必须是2的精确幂。Twiddle-factors W是通过% W = fi_radix2twiddles(N)其中N =长度(X)。％此版本的算法在阶段之前没有伸缩。。％%参见FI_RADIX2FFT_DEMO, FI_M_RADIX2FFT_WITHSCALING。%的参考:Charles Van Loan，快速傅里叶的计算框架。%变换，SIAM，费城，1992，算法1.6.2，第45页。％版权所有The MathWorks, Inc.N =长度(x);T = log2(n);X = fi_bitreverse(X,n);为L = 2^q;r = n/L;L2 = l /2;为K = 0:(r-1)为j = 0:(L2-1) temp = w(L2-1+j+1) * x(k*L+j+L2+1);x(k*L+j+L2+1) = x(k*L+j+1) - temp;x(k*L+j+1) = x(k*L+j+1) + temp;结束结束结束结束

可视化

为了验证您在MATLAB®中正确地实现了算法，运行一个已知信号通过它，并将结果与MATLAB®FFT函数产生的结果进行比较。

如下图所示，误差在MATLAB®内置FFT函数的容忍范围内，验证您已经正确地实现了算法。

Y = fi_m_radix2fft_algorithm1_6_2(x0, w0);fi_fft_demo_plot(真正的(x0), y, y0, Fs,“双数据”，.．.｛“FFT算法1.6.2”，内置的FFT的})；

将函数转换为使用类型表

要从算法中分离数据类型:

创建数据类型定义表。
修改算法代码以使用该表中的数据类型。

这个例子通过创建不同的文件展示了迭代步骤。实际上，您可以对同一个文件进行迭代更改。

原始类型表

使用一个结构创建一个类型表，其中包含变量的原型设置为它们的原始类型。使用基线类型来验证初始转换是否正确，并以编程方式在浮点类型和定点类型之间切换函数。索引变量由MATLAB®Coder™自动转换为整数，因此您不需要在表中指定它们的类型。

将原型值指定为空([])，因为使用的是数据类型，而不是值。

函数T = fi_m_radix2fft_original_types()%FI_M_RADIX2FFT_ORIGINAL_TYPES类型表示例版权所有2015 The MathWorks, Inc.T.x = double([]);T.w = double([]);T.n = double([]);结束

类型感知算法函数

添加类型表T作为函数的输入，并使用它将变量强制转换为特定类型，同时保持算法主体不变。

函数x = fi_m_radix2fft_algorithm1_6_2_typed(x, w, T)%FI_M_RADIX2FFT_ORIGINAL_TYPED基数-2 FFT示例。% Y = FI_M_RADIX2FFT_ALGORITHM1_6_2_TYPED(X, W, T)计算基数-2假设输入向量X的FFT为%复杂。％向量X的长度必须是2的精确幂。Twiddle-factors W是通过% W = fi_radix2twiddles(N)其中N =长度(X)。％% T是一个类型表，用于将变量转换为特定类型，同时保留%的算法体不变。％此版本的算法在阶段之前没有伸缩。。％%参见FI_RADIX2FFT_DEMO, FI_M_RADIX2FFT_WITHSCALING。％%的参考:Charles Van Loan，快速傅里叶的计算框架。%变换，SIAM，费城，1992，算法1.6.2，第45页。版权所有2015 The MathWorks, Inc.％% # codegenN =长度(x);T = log2(n);x = fi_bitreverse_typed(x,n,T);LL = cast(2.^(1:t)，“喜欢”, T.n);rr = cast(n /LL，“喜欢”, T.n);LL2 = cast(LL./2，“喜欢”, T.n);为L = LL(q);R = rr(q);L2 = LL2(q);为K = 0:(r-1)为j = 0:(L2-1) temp = w(L2-1+j+1) * x(k*L+j+L2+1);x(k*L+j+L2+1) = x(k*L+j+1) - temp;x(k*L+j+1) = x(k*L+j+1) + temp;结束结束结束结束

类型感知的位反转函数

添加类型表T作为函数的输入，并使用它将变量强制转换为特定类型，同时保持算法主体不变。

函数x = fi_bitreverse_typed(x,n0,T)%FI_BITREVERSE_TYPED对输入进行反向。% X = FI_BITREVERSE_TYPED(X, n,T)位反转输入序列X，其中% N =长度(X)。％% T是一个类型表，用于将变量转换为特定类型，同时保留%的算法体不变。％%参见FI_RADIX2FFT_DEMO。版权所有The MathWorks, Inc.％% # codegenN = cast(n0，“喜欢”, T.n);Nv2 = bitsra(n,1);J = cast(1，“喜欢”, T.n);为I = 1:(n-1)如果I < j temp = x(j);X (j) = X (i);X (i) = temp;结束K = nv2;而K < j j(:) = j- K;K = bitsra(K,1);结束J (:) = J +k;结束

验证修改后的函数

每次修改函数时，验证结果是否仍然与基线匹配。因为您使用了类型表中的原始类型，所以输出应该是相同的。这验证了您正确地进行了将类型与算法分离的转换。

T1 = fi_m_radix2fft_original_types();获取表中声明的原始数据类型X = cast(x0，“喜欢”, T1.x);W = cast(w0，“喜欢”, T1.w);y = fi_m_radix2fft_algorithm1_6_2_type (x, w, T1);fi_fft_demo_plot(真正的(x), y, y0, Fs,“双数据”，.．.｛“FFT算法1.6.2”，内置的FFT的})；

图中包含3个轴对象。axis对象1包含一个类型为line的对象。该对象表示双数据。坐标轴对象2包含两个line类型的对象。这些对象代表FFT算法1.6.2，内置FFT。Axes对象3包含一个类型为line的对象。该对象表示abs(错误)。

创建一个定点类型表

使用带有变量原型的结构创建定点类型表。将原型值指定为空([])，因为使用的是数据类型，而不是值。

函数T = fi_m_radix2fft_fixed_types()%FI_M_RADIX2FFT_FIXED_TYPES示例函数版权所有2015 The MathWorks, Inc.T.x = fi([]，1,16,14);选择以下类型以确保T.w = fi([]，1,16,14);%的输入具有最大的精度，而不会%溢出T.n = int32([]);%选择int32作为n是一个索引结束

确定定点问题

现在，尝试将输入数据转换为定点，看看算法是否看起来仍然不错。方法使用带符号定点数据的所有默认值fi构造函数。

T2 = fi_m_radix2fft_fixed_types();获取表中声明的定点数据类型X = cast(x0，“喜欢”, T2.x);W = cast(w0，“喜欢”, T2.w);

用定点输入重新运行相同的算法。

y = fi_m_radix2fft_algorithm1_6_2_type (x,w,T2);fi_fft_demo_plot(真正的(x), y, y0, Fs,“定点数据”，.．.｛“定点FFT算法1.6.2”，内置的FFT的})；

注意，定点FFT的幅值图(中心)与内置FFT的图不相似。错误(底部图)比您期望看到的四舍五入错误要大得多，因此很可能发生了溢出。

使用最小/最大仪表来识别溢出

方法从MATLAB®函数中创建一个MEX函数来检测MATLAB®代码buildInstrumentedMex命令。输入buildInstrumentedMex和的输入相同吗fiaccel,但buildInstrumentedMex没有fi对象的限制。的输出buildInstrumentedMex是一个插入插装的MEX函数，因此当MEX函数运行时，将记录所有命名变量和中间值的模拟最小值和最大值。

的“o”option用于为生成的MEX函数命名。如果“o”选项未使用，则MEX函数为MATLAB®函数的名称“_mex”附加。您也可以将MEX函数命名为与MATLAB®函数相同的名称，但您需要记住，MEX函数优先于MATLAB®函数，因此对MATLAB®函数的更改将不会运行，直到重新生成MEX函数，或删除并清除MEX函数。

使用缩放的双数据类型创建输入，这样它的值将达到完整的范围，并且可以识别潜在的溢出。

函数T = fi_m_radix2fft_scaled_fixed_types()%FI_M_RADIX2FFT_SCALED_FIXED_TYPES示例函数版权所有2015 The MathWorks, Inc.DT =“ScaledDouble”；用于fi的数据类型%的构造函数T.x = fi([]，1,16,14，“数据类型”, DT);选择以下类型T.w = fi([]，1,16,14，“数据类型”, DT);保证输入有%最大精度，而不会%溢出T.n = int32([]);%选择int32作为n是一个索引结束

T3 = fi_m_radix2fft_scaled_fixed_types();获取表中声明的定点数据类型X_scaled_double = cast(x0，“喜欢”, T3.x);W_scaled_double = cast(w0，“喜欢”, T3.w);buildInstrumentedMexfi_m_radix2fft_algorithm1_6_2_typed.．.- offt_instrumentedarg游戏{x_scaled_double w_scaled_double T3}

运行仪表MEX函数记录最小/最大值。

y_scaled_double = fft_instrumentation (x_scaled_double,w_scaled_double,T3);

显示检测结果。

showInstrumentationResultsfft_instrumented

您可以从插装结果中看到，在对变量赋值时存在溢出x．

修改算法以解决定点问题

FFT中单个容器的大小最多增长n倍，其中n是FFT的长度。因此，通过将数据扩展为1/n，可以防止任何输入发生溢出。当你只将输入缩放到长度为n的FFT的第一阶段1/n时，你会得到一个与n^2成比例的信噪比[Oppenheim & Schafer 1989，方程9.101]，[Welch 1969]。然而，如果你将FFT的每个阶段的输入缩放1/2，你可以得到1/n的总体缩放，并产生与n成比例的信噪比[Oppenheim & Schafer 1989，方程9.105]，[Welch 1969]。

在定点上缩放1/2的一种有效方法是右移数据。要做到这一点，可以使用位右移算术函数bitsra．在缩放FFT的每个阶段，并优化索引变量计算后，你的算法变成:

函数x = fi_m_radix2fft_withscaling_typed(x, w, T)@ FI_M_RADIX2FFT_WITHSCALING radical -2 FFT示例，每个阶段都有缩放。% Y = FI_M_RADIX2FFT_WITHSCALING_TYPED(X, W, T)计算的基数-2 FFT%输入向量X加上旋转因子W，每个阶段缩放1/2。假设输入X是复数。％向量X的长度必须是2的精确幂。Twiddle-factors W是通过% W = fi_radix2twiddles(N)其中N =长度(X)。％% T是一个类型表，用于将变量转换为特定类型，同时保留%的算法体不变。％此版本的算法在阶段之前没有伸缩。。％%参见FI_RADIX2FFT_DEMO。%的参考:Charles Van Loan，快速傅里叶的计算框架。%变换，SIAM，费城，1992，算法1.6.2，第45页。％版权所有The MathWorks, Inc.％% # codegenN =长度(x);T = log2(n);X = fi_bitreverse(X,n);将索引变量生成为整数常量，因此不计算它们%循环。LL = cast(2.^(1:t)，“喜欢”, T.n);rr = cast(n /LL，“喜欢”, T.n);LL2 = cast(LL./2，“喜欢”, T.n);为L = LL(q);R = rr(q);L2 = LL2(q);为K = 0:(r-1)为j = 0:(L2-1) temp = w(L2-1+j+1) * x(k*L+j+L2+1);x(k*L+j+L2+1) = bitsra(x(k*L+j+1) - temp, 1);x(k*L+j+1) = bitsra(x(k*L+j+1) + temp, 1);结束结束结束结束

用定点数据运行伸缩算法。

X = cast(x0，“喜欢”, T3.x);W = cast(w0，“喜欢”, T3.w);y = fi_m_radix2fft_withscaling_typed(x,w,T3);fi_fft_demo_plot(real(x)， y, y0/n, Fs，“定点数据”，.．.｛“带尺度的定点FFT”，内置的FFT的})；

图中包含3个轴对象。axis对象1包含一个类型为line的对象。该对象表示定点数据。坐标轴对象2包含两个line类型的对象。这些对象代表定点FFT算法1.6.2，内置FFT。Axes对象3包含一个类型为line的对象。该对象表示abs(错误)。

图中包含3个轴对象。axis对象1包含一个类型为line的对象。该对象表示定点数据。坐标轴对象2包含两个line类型的对象。这些对象代表定点FFT与缩放，内置FFT。Axes对象3包含一个类型为line的对象。该对象表示abs(错误)。

您可以看到，缩放的定点FFT算法现在与内置FFT匹配到16位定点数据预期的容错。

清理

运行以下代码恢复全局状态。

fipref (FIPREF_STATE);

参考文献

查尔斯·范洛安，快速傅里叶变换的计算框架，暹罗,1992年。

克里夫硅藻土,MATLAB数值计算，傅里叶分析，2004，第8章。

艾伦·v·奥本海姆和罗纳德·w·谢弗，离散时间信号处理，普伦蒂斯霍尔，1989年。

Peter D. Welch，“定点快速傅立叶变换误差分析”，《IEEE音频与电声汇刊》，AU-17卷，第2期，1969年6月，第151-157页。