stm32g474使用fft
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了stm32g474使用fft相关的知识,希望对你有一定的参考价值。
参考技术A 使用在STM32中对信号进行FFT运算首先需要在程序文件中添加DSP库并将其文件夹加入程序中其中包括头文件与arm_cortexM4lf_math.lib在keil中需要进行如下操作。
使用 Numpy 手动反转 FFT
【中文标题】使用 Numpy 手动反转 FFT【英文标题】:Manually inverting FFT using Numpy 【发布时间】:2016-03-02 14:52:55 【问题描述】:我有一个小脚本用于计算方波的傅立叶变换,当我使用 numpy.fft.ifft()
反转 fft
时,它运行良好并正确返回方波。但是,在将谐波乘以我从numpy.fft.fft()
获得的相应系数之后,我无法通过手动相加谐波来反转变换。下面是我的脚本,我相信你会看到我的意图。
from numpy import zeros, concatenate, sin, pi, linspace
from numpy.fft import fft, fftfreq, ifft
import numpy as np
import matplotlib.pyplot as plt
N = 1024 # samples
T = 1 # period
dt = T/N # sampling period
fs = 1/dt # sampling frequency
t = linspace(0, T, N) # time points
functime = .... # square wave
funcfft = fft(functime) # fft
fftcoeffs = np.abs(funcfft)/N # coefficients, divide by N to get actual coeff.s(I believe?)
freqs = fftfreq(N, dt) # frequencies
plt.plot(freqs, fftcoeffs) # gives me reasonable output
plt.show()
FF = ifft(funcfft)
plt.plot(t, FF) # plots exactly the same function as functime defined above
plt.show()
到目前为止一切都很好。现在我的问题是,如果我在上面的脚本之后运行下面的脚本,我不应该收敛到原始函数吗?:
FFF = zeros(N)
for i in range(300):
FFF += fftcoeffs[i]*sin(2*pi*freqs[i]*t)
plt.plot(t, FFF)
plt.show()
假设range(300)
足以收敛。现在当我这样做时,FFF 与我原来的功能不同。我认为如果我将各个频率的谐波乘以它们对应的系数,我认为这些系数存储在 fftcoeffs 中,然后我会收敛到原始函数。我做错了什么?
更新:根据 DanielSank 的建议,我更新了我的 for 循环如下,不幸的是它没有给我想要的结果:
freqs2 = np.abs(freqs)
freqs2 = np.sort(freqs2)
for k in range(300):
FFF += fftcoeffs[k]*exp(2j*pi*freqs2[k]*t/N)
我不确定我是否在此处执行“按绝对值排序 fftfreq”部分。
【问题讨论】:
您没有正确排序频率 :D 另外,一般来说,当您在任何在线论坛上发布代码时,请包含一个自包含的工作示例 (SCWE)。这被认为是标准做法。不这样做会非常非常频繁地使本来会给您答案的人忽略该帖子。 【参考方案1】:术语
此问题中的任何内容均不针对 fast Fourier transform (FFT)。 FFT 是一种计算discrete Fourier transform (DFT) 的特殊算法,所以我会说“DFT”而不是“FFT”。
在这个答案中,m
表示离散时间索引,k
表示离散频率索引。
什么是 DFT?
这里有几个问题,所有这些问题都源于对 DFT 工作原理的误解。
取自numpy.fft
模块文档字符串,numpy 将离散傅里叶变换定义为
A_k = \sum_m=0^n-1 a_m \exp[-2 \pi i (m k / n)]
LaTeX 表示法表示离散傅里叶变换是复指数的线性组合 exp[2 pi i m k / n]
其中n
是点总数,m
是离散时间索引。
在您的符号中,这将是 exp[2 pi i m k / N]
,因为您使用 N
来表示总点数。
使用exp
,而不是sine
请注意,DFT 使用指数;这些是不是 sine
函数。
如果要从离散傅立叶变换系数构建时域信号,则必须使用 DFT 本身使用的相同函数!
因此,你可以试试这个:
FFF = zeros(N)
for i in range(300):
FFF += fftcoeffs[i]*np.exp(2*pi*i*freqs[i]*t)
plt.plot(t, FFF)
plt.show()
但是,这也会以一种可能让您困惑的方式失败。
别名
最后一块拼图与称为混叠的效果有关。
假设您对信号exp[2 pi i (N + p) m / N]
进行 DFT。
如果你计算它,你会发现所有的A_k
都是零,除了A_p
。
事实上,如果你参加了exp[2 pi i p m / N]
的 DFT,你会得到相同的。
您可以看到频率大于N
的任何复指数都显示为好像它是较低的频率。
特别是,频率为q + b N
的任何复指数,其中b
是任何整数,看起来它的频率为q
。
现在假设我们有一个时域信号cos(2 pi p m / N)
。
这等于
(1/2)[ (exp(2 pi i p m / N) + exp(-2 pi i p m / N) ]
.
负频率对 DFT 产生了有趣的影响。
频率-p
可以写成(N-p) + N
。
这具有q + b N
与q = N - p
和b=1
的形式。
所以,那个负频率-p
最终看起来像N - p
!
numpy 函数fftfreq
知道这一点。
看看fftfreq
的输出,您会看到它从零开始,运行到采样频率的一半(称为奈奎斯特频率),然后变为负数!
这是为了帮助您处理我们刚刚描述的锯齿效应。
所有这一切的结果是,如果您想通过对最低频率傅立叶分量求和来逼近一个函数,您确实不想从fftfreq
中获取最低的几个元素。
相反,您希望按绝对值对 fftfreq 进行排序,然后将复指数与这些频率相加。
还可以查看np.fft.hfft
。
该函数旨在处理实值函数以及与之相关的别名问题。
代码
由于这是一个很难口头讨论的问题,这里有一个脚本可以完全满足您的需求。 请注意,我将 cmets 置于 那些 cmets 描述的代码块之后。 确保您已安装 matplotlib(在您的虚拟环境中......您使用虚拟环境,对吗?)。 如果您有任何问题,请发表评论。
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
pi = np.pi
def square_function(N, square_width):
"""Generate a square signal.
Args:
N (int): Total number of points in the signal.
square_width (int): Number of "high" points.
Returns (ndarray):
A square signal which looks like this:
|____________________
|<-- square_width -->
| ______________
|
|^ ^ ^
index |0 square_width N-1
In other words, the output has [0:N]=1 and [N:]=0.
"""
signal = np.zeros(N)
signal[0:square_width] = 1
return signal
def check_num_coefficients_ok(N, num_coefficients):
"""Make sure we're not trying to add more coefficients than we have."""
limit = None
if N % 2 == 0 and num_coefficients > N // 2:
limit = N/2
elif N % 2 == 1 and num_coefficients > (N - 1)/2:
limit = (N - 1)/2
if limit is not None:
raise ValueError(
"num_coefficients is but should not be larger than ".format(
num_coefficients, limit))
def test(N, square_width, num_coefficients):
"""Test partial (i.e. filtered) Fourier reconstruction of a square signal.
Args:
N (int): Number of time (and frequency) points. We support both even
and odd N.
square_width (int): Number of "high" points in the time domain signal.
This number must be less than or equal to N.
num_coefficients (int): Number of frequencies, in addition to the dc
term, to use in Fourier reconstruction. This is the number of
positive frequencies _and_ the number of negative frequencies.
Therefore, if N is odd, this number cannot be larger than
(N - 1)/2, and if N is even this number cannot be larger than
N/2.
"""
if square_width > N:
raise ValueError("square_width cannot be larger than N")
check_num_coefficients_ok(N, num_coefficients)
time_axis = np.linspace(0, N-1, N)
signal = square_function(N, square_width)
ft = np.fft.fft(signal)
reconstructed_signal = np.zeros(N)
reconstructed_signal += ft[0] * np.ones(N)
# Adding the dc term explicitly makes the looping easier in the next step.
for k in range(num_coefficients):
k += 1 # Bump by one since we already took care of the dc term.
if k == N-k:
reconstructed_signal += ft[k] * np.exp(
1.0j*2 * pi * (k) * time_axis / N)
# This catches the case where N is even and ensures we don't double-
# count the frequency k=N/2.
else:
reconstructed_signal += ft[k] * np.exp(
1.0j*2 * pi * (k) * time_axis / N)
reconstructed_signal += ft[N-k] * np.exp(
1.0j*2 * pi * (N-k) * time_axis / N)
# In this case we're just adding a frequency component and it's
# "partner" at minus the frequency
reconstructed_signal = reconstructed_signal / N
# Normalize by the number of points in the signal. numpy's discete Fourier
# transform convention puts the (1/N) normalization factor in the inverse
# transform, so we have to do it here.
plt.plot(time_axis, signal,
'b.', markersize=20,
label='original')
plt.plot(time_axis, reconstructed_signal.real,
'r-', linewidth=3,
label='reconstructed')
# The imaginary part is zero anyway. We take the real part to
# avoid matplotlib warnings.
plt.grid()
plt.legend(loc='upper right')
【讨论】:
感谢您的回复。首先,我真的不知道如何定义i
又名-1**0.5
。显然-1**0.5
没有定义它。但无论如何,我做了freqs2 = np.abs(freqs)
,然后freqs2 = np.sort(freqs2)
,然后按照你的建议更新了我的for循环,使用exp
而不是sin
。我使用k
作为for 循环变量,而不是i
,这样它就不会与想象中的i
混淆。我猜你忘了在exp()
中包含/N
,我尝试了使用和不使用它,但我在所有尝试中都得到了更糟糕的结果。我怀疑我是否完全理解您的建议。
@Deniz 很难理解代码的书面描述。您可以编辑原始帖子以显示您尝试过的内容吗?
我刚刚看到您对这篇文章的更新并试用了您的脚本。它完美无缺。就我而言,一个周期的正方形(实际上是矩形)信号不像对称的 +a/- 一种正方形信号,所以我必须根据我的情况进行调整,但您的代码仍然可以帮助我理解程序和 DFT 理论。非常感谢您帮助我丹尼尔。我选择这个答案作为正确答案。
@Deniz 很高兴能帮上忙。我发现关于信号处理的文献不能令人满意,所以我就这些主题写了自己的一套笔记。如果您愿意,我很乐意给您一份副本。如果您查看我的物理堆栈交换用户资料并点击我的专业网站的链接,您可以联系我。
@Deniz 你不能@我,因为你在评论我的帖子。用户会自动在自己的帖子中收到有关 cmets 的通知。另外,我想提请您注意,我为您提供了一个完整的工作代码示例,它对您有所帮助。请在您的问题中也这样做。举个例子,问题更容易理解;描述代码几乎是无用的。此外,编写一个最小的示例可以让您在 99% 的时间里解决自己的问题。【参考方案2】:
最大的问题是您只使用了 FFT 结果的幅度 (np.abs),因此丢弃了所有相位信息。
如果您保留 FFT 的复杂相位结果,并在正弦波重新合成中使用该相位信息,您手动添加的谐波将更接近原始结果。
可能的提示:您可能必须在 FFT 之前对波形进行 fftshift 以获得有用的相位结果,具体取决于方波的频率与 FFT 孔径宽度的关系。
【讨论】:
感谢您的回答。使用我上面给出的相同 for 循环,唯一的例外是使用funcfft
而不是 fftcoeffs
,我得到了更糟糕的结果。使用我当前的代码,它给我的完全错误,但它是一个正方形,当我不使用 fftcoeffs
时,它给我的是指数而不是正方形。
为了从复指数中得到一个真实的结果,所有的频率都必须用复共轭来镜像。请注意,完整的复数 FFT 结果也反映了复数共轭。您需要所有这些共轭来获得逆。以上是关于stm32g474使用fft的主要内容,如果未能解决你的问题,请参考以下文章
STM32G431 CubeMx HAL库——使用IDLE+RXNE实现串口不定长数据接收
STM32G070RBT6基于STM32CubeMX创建EXTI外部中断工程
STM32G070RBT6基于Arduino框架下eeprom使用示例