频域信号处理

最新推荐文章于 2022-12-15 11:58:40 发布

sinat_18131557

最新推荐文章于 2022-12-15 11:58:40 发布

阅读量2.6k

点赞数 2

分类专栏：频域

本文链接：https://blog.csdn.net/sinat_18131557/article/details/103431977

版权

频域专栏收录该内容

1 篇文章 1 订阅

订阅专栏

频域信号处理

FFT

FFT变换是针对一组数值进行运算的，这组数的长度N必须是2的整数次幂，例如64, 128, 256等等；数值可以是实数也可以是复数，通常我们的时域信号都是实数，因此下面都以实数为例。我们可以把这一组实数想像成对某个连续信号按照一定取样周期进行取样而得来，如果对这组N个实数值进行FFT变换，将得到一个有N个复数的数组，我们称此复数数组为频域信号，此复数数组符合如下规律：

下标为0和N/2的两个复数的虚数部分为0，
下标为i和N-i的两个复数共轭，也就是其虚数部分数值相同、符号相反。

FFT可以将信号转化为频谱。转换为频域信号之后我们可以很方便地分析出信号的频率成分，在频域上进行处理，最终还可以将处理完毕的频域信号通过IFFT(逆变换)转换为时域信号，实现许多在时域无法完成的信号处理算法。
FFT能返回采样频率一半以下的频率分布情况。原理不过多去介绍如何实现的，介绍以下MATLAB以及Python是怎么实现的。

%matlab
function [ output_args ] = mfft( data,Fs)
    %MFFT 此处显示有关此函数的摘要
    %   此处显示详细说明
    data = data-mean(data); %去除均值，避免频谱图在0点的幅值非常大。
    h = figure;
    subplot(211);
    plot(data);
    subplot(212);
    N = length(data);
    n = 0:N-1;
    y=fft(data,N);
    mag = abs(y);
    f=n*Fs/N;
    plot(f(1:fix(N/2)),mag(1:fix(N/2)),'r');
    xlabel('频率/Hz')
end

这里用matlab的内置函数fft对信号进行变换，并且对坐标进行变化，就能得到变换结果，并且绘图出来。

# python
import numpy as np
from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt
import seaborn

# 采样点选择1400个，因为设置的信号频率分量最高为600Hz，根据采样定理知采样频率要大于信号频率2倍，所以这里设置采样频率为1400Hz（即一秒内有1400个采样点）
x = np.linspace(0, 1, 1400)

# 设置需要采样的信号，频率分量有180，390和600
y = 7 * np.sin(2 * np.pi * 180 * x) + 1.5 * np.sin(2 * np.pi * 390 * x) + 5.1 * np.sin(2 * np.pi * 600 * x)

yy = fft(y)  # 快速傅里叶变换
yreal = yy.real  # 获取实数部分
yimag = yy.imag  # 获取虚数部分

yf = abs(fft(y))  # 取模
yf1 = abs(fft(y)) / (len(x) / 2)  # 归一化处理
yf2 = yf1[range(int(len(x) / 2))]  # 由于对称性，只取一半区间

xf = np.arange(len(y))  # 频率
xf1 = xf
xf2 = xf[range(int(len(x) / 2))]  # 取一半区间

# 原始波形
plt.subplot(221)
plt.plot(x[0:50], y[0:50])
plt.title('Original wave')
# 混合波的FFT（双边频率范围）
plt.subplot(222)
plt.plot(xf, yf, 'r')  # 显示原始信号的FFT模值
plt.title('FFT of Mixed wave(two sides frequency range)', fontsize=7, color='#7A378B')  # 注意这里的颜色可以查询颜色代码表
# 混合波的FFT（归一化）
plt.subplot(223)
plt.plot(xf1, yf1, 'g')
plt.title('FFT of Mixed wave(normalization)', fontsize=9, color='r')

plt.subplot(224)
plt.plot(xf2, yf2, 'b')
plt.title('FFT of Mixed wave)', fontsize=10, color='#F08080')

plt.show()

python是通过scipy.fftpack包实现的，在绘图的224里面，实际就是需要的频谱图。

快速卷积

我们知道，信号x经过系统h之后的输出y是x和h的卷积，虽然卷积的计算方法很简单，但是当x和h都很长的时候，卷积计算是非常耗费时间的。因此对于比较长的系统h，需要找到比直接计算卷积更快的方法。

信号系统理论中有这样一个规律：时域的卷积等于频域的乘积，因此要计算时域的卷积，可以将时域信号转换为频域信号，进行乘积运算之后再将结果转换为时域信号，实现快速卷积。

由于FFT运算可以高效地将时域信号转换为频域信号，其运算的复杂度为 O(N*log(N))，因此三次FFT运算加一次乘积运算的总复杂度仍然为 O(N*log(N)) 级别，而卷积运算的复杂度为 O(N*N)，显然通过FFT计算卷积要比直接计算快速得多。这里假设需要卷积的两个信号的长度都为N。

但是有一个问题：FFT运算假设其所计算的信号为周期信号，因此通过上述方法计算出的结果实际上是两个信号的循环卷积，而不是线性卷积。为了用FFT计算线性卷积，需要对信号进行补零扩展，使得其长度长于线性卷积结果的长度。

例如，如果我们要计算数组a和b的卷积，a和b的长度都为128，那么它们的卷积结果的长度为 len(a) + len(b) - 1 = 257。为了用FFT能够计算其线性卷积，需要将a和b都扩展到256。下面的程序演示这个计算过程：

# -*- coding: utf-8 -*-
import numpy as np

def fft_convolve(a,b):
    n = len(a)+len(b)-1
    N = 2**(int(np.log2(n))+1)
    A = np.fft.fft(a, N)
    B = np.fft.fft(b, N)
    return np.fft.ifft(A*B)[:n]

if __name__ == "__main__":
    a = np.random.rand(128)
    b = np.random.rand(128)
    c = np.convolve(a,b)
    
    print np.sum(np.abs(c - fft_convolve(a,b)))

可以看到误差是非常小的。
在这段程序中，a,b的长度为128，其卷积结果c的长度为n=255，我们通过下面的算式找到大于n的最小的2的整数次幂：

N = 2**(int(np.log2(n))+1)

分段运算

现在考虑对于输入信号x和系统响应h的卷积运算，通常x是非常长的，例如要对某段录音进行滤波处理，假设取样频率为8kHz，录音长度为1分钟的话，那么x的长度为480000。而且x的长度也可能不是固定的，例如我们可能需要对麦克风的连续输入信号进行滤波处理。而h的长度通常都是固定的，例如它是某个房间的冲击响应，或者是某种FIR滤波器。

根据前面的介绍，为了有效地利用FFT计算卷积，我们希望它的两个输入长度相当，于是就需要对信号x进行分段处理。对卷积的分段运算被称作：overlap-add运算。

在这里插入图片描述

原始信号x长度为300，将它分为三段，分别与滤波器系数h进行卷积计算，h的长度为101，因此每段输出200个数据，图中用绿色标出每段输出的200个数据。这3段数据按照时间顺序进行求和之后得到结果和原始信号的卷积是相同的。

因此将持续的输入信号x和滤波器h进行卷积的运算可以按照如下步骤进行，假设h的长度为M：

建立一个缓存，其大小为N+M-1，初始值为0
每次从x中读取N个数据，和h进行卷积，得到N+M-1个数据，和缓存中的数据进行求和，并放进缓存中，然后输出缓存前N个数据
将缓存中的数据向左移动N个元素，也就是让缓存中的第N个元素成为第0个元素，后面的N个元素全部设置为0
跳转到2重复运行

下面是实现这一算法的演示程序：

# -*- coding: utf-8 -*-
import numpy as np

x = np.random.rand(1000)
h = np.random.rand(101)
y = np.convolve(x, h)

N = 50  # 分段大小
M = len(h)  # 滤波器长度

output = []

# 缓存初始化为0
buffer = np.zeros(M + N - 1, dtype=np.float64)

for i in range(len(x) // N):
    # 从输入信号中读取N个数据
    xslice = x[i * N:(i + 1) * N]
    # 计算卷积
    yslice = np.convolve(xslice, h)
    # 将卷积的结果加入到缓冲中
    buffer += yslice
    # 输出缓存中的前N个数据，注意使用copy，否则输出的是buffer的一个视图
    output.append(buffer[:N].copy())
    # 缓存中的数据左移动N个元素
    buffer[0:M - 1] = buffer[N:]
    # 后面的补0
    buffer[M - 1:] = 0

# 将输出的数据组合为数组
y2 = np.hstack(output)
# 计算和直接卷积的结果之间的误差
print(np.sum(np.abs(y2 - y[:len(x)])))

注意第23行需要输出缓存前N个数据的拷贝，否则输出的是数组的一个视图，当此后buffer更新时，视图中的数据会一起更新。

将FFT快速卷积和overlap-add相结合，可以制作出一些快速的实时数据滤波算法。但是由于FFT卷积对于两个长度相当的数组时最为有效，因此在分段时也会有所限制：例如如果滤波器的长度为2048，那么理想的分段长度也为2048，如果将分段长度设置得过低，反而会增加运算量。因此在实时性要求很强的系统中，只能采用直接卷积。

Hilbert变换

Hilbert变换能在振幅保持不变的情况下将输入信号的相角偏移90度，简单地说就是能将正弦波形转换为余弦波形，下面的程序验证这一特性：

# -*- coding: utf-8 -*-
from scipy import fftpack
import numpy as np
import matplotlib.pyplot as pl

# 产生1024点4个周期的正弦波
t = np.linspace(0, 8*np.pi, 1024, endpoint=False)
x = np.sin(t)

# 进行Hilbert变换
y = fftpack.hilbert(x)
pl.plot(x, label=u"原始波形")
pl.plot(y, label=u"Hilbert转换后的波形")
pl.legend()
pl.show()

关于程序的几点说明：

hilbert转换函数在scipy.fftpack函数库中
为了生成完美的正弦波，需要计算整数个周期，因此调用linspace时指定endpoint=False，这样就不会包括区间的结束点：8*np.pi

Hilbert正变换的相角偏移符号:本文中将相角偏移+90度成为Hilbert正变换。有的文献书籍正好将定义倒转过来：将偏移+90度成为Hilbert负变换，而偏移-90度成为Hilbert正变换。

相角偏移90度相当于复数平面上的点与虚数单位1j相乘，因此Hilbert变换的频率响应可以用如下公式给出：
$H(\omega)=j\cdot sgn(\omega)$
其中 $\omega$ 为圆频率，sgn函数为符号函数：
$sgn(\omega)=\left\{\begin{array}{lc}1,&\omega>0\\0,&\omega=0\\-1,&\omega<0\end{array}\right.$

我们可以将其频率响应理解为：

直流分量为0
正频率成分偏移+90度
负频率成分偏移-90度

由于对于实数信号来说，正负频率成分共轭，因此对实数信号进行Hilbert变换之后仍然是实数信号。下面的程序验证Hilbert变换的频率响应：

Hilbert变换可以用于求信号的包络线，实现方式为: $envelope=\sqrt{H(x)^2+x^2}$ 其中x为原始载波波形，H(x)为x的Hilbert变换之后的波形，envelope为信号x的包络。其原理很容易理解：假设x为正弦波，那么H(x)为余弦波，根据公式：
$sin(x)^2+cos(x)^2=1$ 可知envelope恒等于1，为sin(t)信号的包络。下面的程序验证这一算法：

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl
from scipy import fftpack

t = np.arange(0, 0.3, 1/2000.0)
x = np.sin(2*np.pi*100*t) * (np.sin(2*np.pi*10*t) + np.sin(2*np.pi*7*t) + 3.0)
hx = fftpack.hilbert(x)

pl.plot(x, label=u"sig")
pl.plot(np.sqrt(x**2 + hx**2), "r", linewidth=2, label=u"envelope")
pl.title(u"sig & envelope")
pl.legend()
pl.show()

在这里插入图片描述

Hilbert变换在信号的高频和低频处包络计算出现较大的误差。而在中频部分能很好地计算出包络的形状。

参考：

https://blog.csdn.net/qq_39516859/article/details/79766697
https://docs.huihoo.com/scipy/scipy-zh-cn/frequency_process.html

sinat_18131557

关注

2
点赞
踩
12

收藏

觉得还不错? 一键收藏
0
评论
频域信号处理

频域信号处理FFTFFT变换是针对一组数值进行运算的，这组数的长度N必须是2的整数次幂，例如64, 128, 256等等；数值可以是实数也可以是复数，通常我们的时域信号都是实数，因此下面都以实数为例。我们可以把这一组实数想像成对某个连续信号按照一定取样周期进行取样而得来，如果对这组N个实数值进行FFT变换，将得到一个有N个复数的数组，我们称此复数数组为频域信号，此复数数组符合如下规律：下...
复制链接

扫一扫