第七章离散傅里叶变换

最新推荐文章于 2022-08-03 15:06:15 发布

weixin_45965693

最新推荐文章于 2022-08-03 15:06:15 发布

阅读量536

点赞数

分类专栏： Python 数字信号处理

本文链接：https://blog.csdn.net/weixin_45965693/article/details/113867463

版权

Python 同时被 2 个专栏收录

54 篇文章 3 订阅

订阅专栏

数字信号处理

14 篇文章 6 订阅

订阅专栏

离散傅里叶变换（DFT）与离散余弦变换（DCT）之间唯一的区别在于，不再使用余弦函数，而是使用负指数函数。

7.1 复指数

幂可以表示为乘幂序列，这个定义适合实数、虚数
$e^{φ}=1+φ+φ^{2}/2!+φ^{3}/3!+…\\ e^{iφ}=1+iφ-φ^{2}/2!-iφ^{3}/3!+…\\ e^{iφ}=cosφ+i sinφ$
可以理解为 $e^{iφ}$ 是一个大小为1的复数，在复平面上，它总是在单位圆上，也是一个向量，与x轴正半轴的弧度角便是 $φ$
$e^{\alpha+iφ}=e^{\alpha}e^{iφ}=Ae^{iφ}$
其中A是实数，意为振幅， $e^{iφ}$ 是单位复数，意为角度

NumPy提供了一个版本的exp支持复数

phi = 1.5
z = np.exp(1j * phi)#大小为1角度为φ
z.real
z.imag
abs(z)#1.0
np.absolute(z)#1.0
np.angle#1.5

7.2 复信号

$e^{iφ(t)}=cosφ(t)+i sinφ(t)$
复正弦信号 $φ(t)=2\pi ft$
$e^{i2\pi ft}=cos2\pi ft+i sin2\pi ft$
更一般信号可能从相位 $φ_{0}$ 开始 $e^{i(2\pi ft+φ_{0})}$
实现该采样结果为复数的信号

 class ComplexSinusoid(Sinusoid):
	def evaluate(self, ts):
		phases = PI2 * self.freq * ts + self.offset
		ys = self.amp *np.exp(1j * phases)
		return ys

复信号是一序列复数，如何解释它呢？这里有一个不太令人满意的答案：
它包含两个信号——实数部分和虚数部分

7.3 合成问题

给定各个复数元素的频率和振幅，如何对信号求值
方法一：创建一个ComplexSinusoid对象然后叠加它们

def synthesize1(amp, fs, ts):
	components = [thinkdsp.ComplexSinusoid(freq, amp) for amp, freq in zip(amps, fs)]
	signal = thinkdsp.SumSignal(*components)
	ys = signal.evaluate(ts)
	return ys

解释复信号

amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = [100, 200, 300, 400]
framerate = 11025
ts = np.linspace(0, 1, framerate)
ys = synthesize1(amps, fs, ts)
n = 500
thinkplot.plot(ts[:n], ys[:n].real, label='real')
thinkplot.plot(ts[:n], ys[:n].imag, label='imag')

在这里插入图片描述
它们包含的频率元素的对应比例是一样的，对我们耳朵来说，感知是一样的，虽然实数部分的余弦波的和，虚数部分是正弦波的和。

7.4 使用矩阵合成

PI2 = np.pi * 2
def synthesize2(amps, fs, ts):
	args = np.outer(ts, fs)
	M = np.exp(1j * PI2 * args)
	ys = np.dot(M, amps)
	return ys

若复振幅和复正弦相乘
$Ae^{iφ_{0}}·e^{i2\pi ft}=Ae^{i(2\pi ft+φ_{0})}$
等于将振幅乘以A，将相位加上了 $φ_{0}$

不同相位差的波形：

phi = 1.5
amps2 = amps * np.exp(1j * phi)
ys2 = synthesize2(amps2, fs, ts)
thinkplot.plot(ts[:n], ys.real[:n])
thinkplot.plot(ts[:n], ys2.real[:n])

在这里插入图片描述

$φ_{0}=1.5$ 时各个频率元素被移动了大约十分之一个周期。但是不同频率元素的周期不同，把这些元素加和时，波形就看其了不同了

因为复振幅影响相位，相位影响波形，所以处理复振幅问题，我们就可以解决合成问题

7.5 分析问题

分析问题就是合成问题的逆命题
方法一：解线性方程组 $M a = y$

def analyze1(ys, fs, ts):
	args = np.outer(ts, fs)
	M = np.exp(1j * PI2 * args)
	amps = np.linalg.solve(M,ys)
	return amps

7.6 快速分析

DCT提速是通过选择fs和ts让M正交
DFT也是样的思路，不过对于M复数矩阵，我们需要它的酉矩阵

M的逆是M的共轭转置（转置并去每个元素虚部的负值）

N = 4
ts = np.arange(N) / N
fs = np.arange(N)
args = np.outer(ts, fs)
M = np.exp(1j * PI2 * args)
#测试M是否为酉矩阵
MstarM = M.conj().transpose().dot(M)#4I有一个额外因子N

def analyze2(ys, fs, ts):
	args = np.outer(ts, fs)
	M = np.exp(1j * PI2 * args)
	amps = M.conj().transpose().dot(ys) / N
	return amps

7.7 DFT

重写analyze，让它只接受ys，从而自行算出freq和fs

#创造一个函数来计算合成矩阵M
def synthesis_matrix(N):
	ts = np.arange(N) / N
	fs = np.arange(N)
	args = np.outer(ts, fs)
	M = np.exp(1j * PI2 * args)
	return M
def analyze3(ys):	
	N = len(ys)
	M = sythesis_matrix(N)
	amps = M.conj().transpose().dot(ys) / N
	return amps
def dft(ys):
	N = len(ys)
	M = synthesis_matrix(N)	
	amps = M.conj().transpose().dot(ys)
	return amps
def idft(ys):
	N = len(ys)
	M = synthesis_matrix(N)
	amps = M.dot(ys) / N 
	return amps

7.8 DFT是周期性的

$\sum_{n}y[n]exp(-2\pi ink / N)\\ DFT(y)[k+N] = \sum_{n}y[n]exp(-2\pi in(k+N) / N)\\ DFT(y)[k+N] = \sum_{n}y[n]exp(-2\pi ink / N)exp(-2\pi inN / N)\\ DFT(y)[k+N] = \sum_{n}y[n]exp(-2\pi ink / N)$

DFT的每一个元素是一个相关，它度量的是波形数组和一个复指数在特定频率上的相似性。

7.9 实信号

np.ftt.rftt只支持实信号
本章的DFT更一般化，把这个“完全版DFT”应用到实信号：

signal = thinkdsp.SawtoothSignal(freq=500)
wave = signal.make_wave(duration=0.1,framerate=10000)
hs = dft(wave.ys)#包含波形的复DFT
amps = np.absolute(hs)#包含各个频率的振幅
#但是这些振幅对应的频率是什么呢?
#dft不知道采样率，它假设波形的持续时间为一个时间单位，假设采样率为每单位时间N次
fs = np.arange(N)
#为了得到频率，需要将随意的时间单位转换为秒
#通过这个转化，频率的范围变成了从0到实际的帧率10kHz
fs = np.arange(N) * framerate / N
#绘制频谱
thinkplot.plot(fs, amps)
thinkplot.config(xlabel='frequency (Hz)',
		 ylabel='amplitude')