信号生成及DFT的python实现方式

DFT

DFT(Discrete Fourier Transform),离散傅里叶变化,可以将离散信号变换到频域,它的公式非常简单:

离散频率下标为k时的频率大小

离散时域信号序列

信号序列的长度,也就是采样的个数

如果你刚接触DFT,并且之前没有信号处理的相关经验,那么第一次看到这个公式,你可能有一些疑惑,为什么这个公式就能进行时域与频域之间的转换呢?

这里,我不打算去解释它,因为我水平有限,说的不清楚。相反,在这里我想介绍,作为一个程序员,如何如实现DFT

从矩阵的角度看DFT

DFT的公式,虽然简单,但是理解起来比较麻烦,我发现如果用矩阵相乘的角度来理解上面的公式,就会非常简单,直接上矩阵:

OK,通过上面的表示,我们很容易将DFT理解成为一种矩阵相乘的操作,这对于我们编码是很容易的。

Talk is cheap, show me the code

根据上面的理解,我们只需要构建出S SS矩阵,然后做矩阵相乘,就等得到DFT的结果

在这之前,我们先介绍如何生成正弦信号,以及如何用scipy中的fft模块进行DFT操作,以验证我们的结果是否正确

正弦信号

A: 幅度

f: 信号频率

n: 时间下标

T: 采样间隔, 等于 1/fs,fs为采样频率

ϕ \phiϕ: 相位

下面介绍如何生成正弦信号

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
def generate_sinusoid(N, A, f0, fs, phi):
 '''
 N(int) : number of samples
 A(float) : amplitude
 f0(float): frequency in Hz
 fs(float): sample rate
 phi(float): initial phase

 return
 x (numpy array): sinusoid signal which lenght is M
 '''

 T = 1/fs
 n = np.arange(N) # [0,1,..., N-1]
 x = A * np.cos( 2*f0*np.pi*n*T + phi )

 return x

N = 511
A = 0.8
f0 = 440
fs = 44100
phi = 0

x = generate_sinusoid(N, A, f0, fs, phi)

plt.plot(x)
plt.show()

# 另一种生成正弦信号的方法,生成时长为t的序列
def generate_sinusoid_2(t, A, f0, fs, phi):
 '''
 t (float) : 生成序列的时长
 A (float) : amplitude
 f0 (float) : frequency
 fs (float) : sample rate
 phi(float) : initial phase

 returns
 x (numpy array): sinusoid signal sequence
 '''

 T = 1.0/fs
 N = t / T

 return generate_sinusoid(N, A, f0, fs, phi)

A = 1.0
f0 = 440
fs = 44100
phi = 0
t = 0.02

x = generate_sinusoid_2(t, A, f0, fs, phi)

n = np.arange(0, 0.02, 1/fs)
plt.plot(n, x)

Scipy FFT

介绍如何Scipy的FFT模块计算DFT

注意,理论上输入信号的长度必须是才能做FFT,而scipy中FFT却没有这样的限制

这是因为当长度不等于时,scipy fft默认做DFT

from scipy.fftpack import fft

# generate sinusoid
N = 511
A = 0.8
f0 = 440
fs = 44100
phi = 1.0
x = generate_sinusoid(N, A, f0, fs, phi)

# fft is
X = fft(x)
mX = np.abs(X) # magnitude
pX = np.angle(X) # phase

# plot the magnitude and phase
plt.subplot(2,1,1)
plt.plot(mX)

plt.subplot(2,1,2)
plt.plot(pX)
plt.show()

自己实现DFT

自己实现DFT的关键就是构造出S,有两种方式:

def generate_complex_sinusoid(k, N):
 '''
 k (int): frequency index
 N (int): length of complex sinusoid in samples

 returns
 c_sin (numpy array): the generated complex sinusoid (length N)
 '''

 n = np.arange(N)

 c_sin = np.exp(1j * 2 * np.pi * k * n / N)

 return np.conjugate(c_sin)

def generate_complex_sinusoid_matrix(N):
 '''
 N (int): length of complex sinusoid in samples

 returns
 c_sin_matrix (numpy array): the generated complex sinusoid (length N)
 '''

 n = np.arange(N)
 n = np.expand_dims(n, axis=1)  # 扩充维度,将1D向量,转为2D矩阵,方便后面的矩阵相乘

 k = n

 m = n.T * k / N     # [N,1] * [1, N] = [N,N]

 S = np.exp(1j * 2 * np.pi * m)  # 计算矩阵 S

 return np.conjugate(S)
# 生成信号,用于测试
N = 511
A = 0.8
f0 = 440
fs = 44100
phi = 1.0
x = generate_sinusoid(N, A, f0, fs, phi)

# 第一种方式计算DFT
X_1 = np.array([])
for k in range(N):
 s = generate_complex_sinusoid(k, N)
 X_1 = np.append(X_1, np.sum(x * s))

mX = np.abs(X_1)
pX = np.angle(X_1)

# plot the magnitude and phase
plt.subplot(2,1,1)
plt.plot(mX)

plt.subplot(2,1,2)
plt.plot(pX)
plt.show()

# 结果和scipy的结果基本相同

# 第二种方法计算DFT
S = generate_complex_sinusoid_matrix(N)
X_2 = np.dot(S, x)

mX = np.abs(X_2)
pX = np.angle(X_2)

# plot the magnitude and phase
plt.subplot(2,1,1)
plt.plot(mX)

plt.subplot(2,1,2)
plt.plot(pX)
plt.show()

总结

回顾了DFT的计算公式,并尝试用矩阵相乘的角度来理解DFT

介绍了两种生成正弦信号的方法

实现了两种DFT的计算方法

完整代码在这里

以上这篇信号生成及DFT的python实现方式就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持我们。

(0)

相关推荐

  • python实现随机加减法生成器

    本文实例为大家分享了python实现随机加减法生成器的具体代码,供大家参考,具体内容如下 为了让外甥女练习算术,用python给她写了个自动出加减法的小程序. 该程序使用了文字转语音的库pyttsx,程序运行时,会有相对应的语音提示.pyttsx文档 为了防止小孩乱按键盘,导致非法输入,我添加了异常处理: def validate(num):#判断输入是否非法 try: num=int(num) except: say("非法输入,请重新输入") return False return

  • python 实现快速生成连续、随机字母列表

    0.摘要 本文介绍了生成连续和随机字母表的方法,用于快速生成大量字母数据. 主要使用chr()函数,将数字通过ASCII表转换为相应字母. 1.chr() 函数 chr() 用一个范围在 range(256)内的(就是0-255)整数作参数,返回一个对应的字符. 输入:可以是10进制也可以是16进制的形式的数字. print(chr(48), chr(49), chr(97)) # 十进制 #result:0 1 a print(chr(0x30), chr(0x31), chr(0x61))

  • python随机生成库faker库api实例详解

    废话不多说,直接上代码! # -*- coding: utf-8 -*- # @Author : FELIX # @Date : 2018/6/30 9:49 from faker import Factory # zh_CN 表示中国大陆版 fake = Factory().create('zh_CN') # 产生随机手机号 print(fake.phone_number()) # 产生随机姓名 print(fake.name()) # 产生随机地址 print(fake.address())

  • python生成任意频率正弦波方式

    如下所示: def signal_xHz(A, fi, time_s, sample): return A * np.sin(np.linspace(0, fi * time_s * 2 * np.pi , sample* time_s)) A:为信号幅值 fi:为信号频率 time_s:为时间长度(s) sample:为信号采样频率 补充拓展:Python FFT合成波形实例 使用Python numpy模块带的FFT函数合成矩形波和方波,增加对离散傅里叶变换的理解. 导入模块 import

  • 信号生成及DFT的python实现方式

    DFT DFT(Discrete Fourier Transform),离散傅里叶变化,可以将离散信号变换到频域,它的公式非常简单: 离散频率下标为k时的频率大小 离散时域信号序列 信号序列的长度,也就是采样的个数 如果你刚接触DFT,并且之前没有信号处理的相关经验,那么第一次看到这个公式,你可能有一些疑惑,为什么这个公式就能进行时域与频域之间的转换呢? 这里,我不打算去解释它,因为我水平有限,说的不清楚.相反,在这里我想介绍,作为一个程序员,如何如实现DFT 从矩阵的角度看DFT DFT的公式

  • python生成并处理uuid的实现方式

    UUID(Universally Unique Identifier)是通用唯一识别码,在许多领域用作标识,比如我们常用的数据库也可以用它来作为主键,原理上它是可以对任何东西进行唯一的编码的. 作为新手一看到类似varchar(40)这样的主键就觉得有点蒙圈了,字符串型也不能自增啊,这里就应该应用UUID了. 数据库一般都有自己的办法生成UUID,但虽然可以用,但这玩意考虑到可读性和有点坑的长度还是尽量不要用这玩意做主键···咳,有点跑题··· 下面就简单说明一下python是如何生成UUID的

  • 使用Python项目生成所有依赖包的清单方式

    1.安装所需工具 pip install pipreqs 2.进入到python项目主目录 pipreqs ./ 3.完成上面命令会生成requirements.txt 4.sudo pip install -r requirements.txt即可 补充知识:解决Python开发过程中依赖库打包问题的方法 在Python开发的过程中,经常会遇到各种各样的小问题,比如在一台计算机上调试好的程序,迁移到另外一台机子上后往往会应为工程项目依赖库的缺失而造成错误. 除了一遍又一遍对着被抛出错误去重新i

  • Python进程间通信方式

    目录 一.通信方式 二.Queue介绍 三.方法介绍 三.生产者和消费者模型 四.什么是生产者消费者模式 实现方式一:Queue 实现方式二:利用JoinableQueue 一.通信方式 进程彼此之间互相隔离,要实现进程间通信(IPC),multiprocessing模块主要通过队列方式 队列:队列类似于一条管道,元素先进先出 需要注意的一点是:队列都是在内存中操作,进程退出,队列清空,另外,队列也是一个阻塞的形态 二.Queue介绍 创建队列的类(底层就是以管道和锁定的方式实现): Queue

  • 基于YUV 数据格式详解及python实现方式

    YUV 数据格式概览 YUV 的原理是把亮度与色度分离,使用 Y.U.V 分别表示亮度,以及蓝色通道与亮度的差值和红色通道与亮度的差值.其中 Y 信号分量除了表示亮度 (luma) 信号外,还含有较多的绿色通道量,单纯的 Y 分量可以显示出完整的黑白图像.U.V 分量分别表示蓝 (blue).红 (red) 分量信号,它们只含有色彩 (chrominance/color) 信息,所以 YUV 也称为 YCbCr,C 意思可以理解为 (component 或者 color). 维基百科上的 RGB

  • K最近邻算法(KNN)---sklearn+python实现方式

    k-近邻算法概述 简单地说,k近邻算法采用测量不同特征值之间的距离方法进行分类. k-近邻算法 优点:精度高.对异常值不敏感.无数据输入假定. 缺点:计算复杂度高.空间复杂度高. 适用数据范围:数值型和标称型. k-近邻算法(kNN),它的工作原理是:存在一个样本数据集合,也称作训练样本集,并且样本集中每个数据都存在标签,即我们知道样本集中每一数据与所属分类的对应关系.输入没有标签的新数据后,将新数据的每个特征与样本集中数据对应的特征进行比较,然后算法提取样本集中特征最相似数据(最近邻)的分类标

  • Mysql数据库反向生成Django里面的models指令方式

    python manage.py inspectdb 或 python manage.py inspect > app/models.py 补充知识:Django框架MySQL数据库到models模型的映射关系 一.前言 我的数据库已经用MySQL Workbench设计好了,也插入了一些测试数据,现在开始在Django中设计models模型.本以为顺风顺水,没想到也遇到一些bug,现在记录一下踩坑填坑过程. 二.设计models模型 1. 如果数据库中表的数量比较多,可以先导出,然后查看对应表

  • python多线程方式执行多个bat代码

    python多线程方式执行多个bat的代码,感兴趣的朋友可以参考下. import threading from win32api import * class MyThread(threading.Thread): def __init__(self, bat_path, **kwargs): threading.Thread.__init__(self, **kwargs) self.bat_path = bat_path def run(self): ShellExecute(0, Non

  • 在AOP中Spring生成代理类的两种方式

    Java 动态代理.具体有如下四步骤: 通过实现 InvocationHandler 接口创建自己的调用处理器: 通过为 Proxy 类指定 ClassLoader 对象和一组 interface 来创建动态代理类: 通过反射机制获得动态代理类的构造函数,其唯一参数类型是调用处理器接口类型: 通过构造函数创建动态代理类实例,构造时调用处理器对象作为参数被传入. 在AOP中,Spring通过生成代理类,来完成切面的织入. Spring生成代理类有2种方式. 如果目标对象实现的是一个接口,Sprin

  • Python多进程方式抓取基金网站内容的方法分析

    本文实例讲述了Python多进程方式抓取基金网站内容的方法.分享给大家供大家参考,具体如下: 在前面这篇//www.jb51.net/article/162418.htm我们已经简单了解了"python的多进程",现在我们需要把抓取基金网站(28页)内容写成多进程的方式. 因为进程也不是越多越好,我们计划分3个进程执行.意思就是 :把总共要抓取的28页分成三部分. 怎么分呢? # 初始range r = range(1,29) # 步长 step = 10 myList = [r[x:

随机推荐