我应该如何将 scipy.fftpack 输出向量相乘?

2023-12-28

The scipy.fftpack.rfft函数将 DFT 作为浮点数向量返回,在实数部分和复数部分之间交替。这意味着要与 DFT 相乘(对于卷积),我将必须“手动”进行复杂的乘法,这看起来相当棘手。这一定是人们经常做的事情 - 我认为/希望有一个我没有发现的简单技巧可以有效地做到这一点?

基本上我想修复此代码,以便两种方法给出相同的答案:

import numpy as np
import scipy.fftpack as sfft

X = np.random.normal(size = 2000)
Y = np.random.normal(size = 2000)
NZ = np.fft.irfft(np.fft.rfft(Y) * np.fft.rfft(X))
SZ = sfft.irfft(sfft.rfft(Y) * sfft.rfft(X))    # This multiplication is wrong

NZ
array([-43.23961083,  53.62608086,  17.92013729, ..., -16.57605207,
     8.19605764,   5.23929023])
SZ
array([-19.90115323,  16.98680347,  -8.16608202, ..., -47.01643274,
    -3.50572376,  58.1961597 ])

注意:我知道 fftpack 包含convolve函数,但我只需要 fft 变换的一半 - 我的过滤器可以提前 fft 一次,然后一遍又一遍地使用。


You don't必须翻转回np.float64 and hstack。您可以创建一个空的目标数组,其形状与sfft.rfft(Y) and sfft.rfft(X),然后创建一个np.complex128视图并用乘法的结果填充该视图。这将根据需要自动填充目标数组。
如果我重新举你的例子:

import numpy as np
import scipy.fftpack as sfft

X = np.random.normal(size = 2000)
Y = np.random.normal(size = 2000)
Xf = np.fft.rfft(X)
Xf_cpx = Xf[1:-1].view(np.complex128)
Yf = np.fft.rfft(Y)
Yf_cpx = Yf[1:-1].view(np.complex128)

Zf = np.empty(X.shape)
Zf_cpx = Zf[1:-1].view(np.complex128)

Zf[0] = Xf[0]*Yf[0]

# the [...] is important to use the view as a reference to Zf and not overwrite it
Zf_cpx[...] = Xf_cpx * Yf_cpx 

Zf[-1] = Xf[-1]*Yf[-1]

Z = sfft.irfft.irfft(Zf)

就是这样! 如果您希望代码更通用并处理奇数长度(如 Jaime 的答案中所述),您可以使用简单的 if 语句。 这是一个可以完成您想要的功能的函数:

def rfft_mult(a,b):
    """Multiplies two outputs of scipy.fftpack.rfft"""
    assert a.shape == b.shape
    c = np.empty( a.shape )
    c[...,0] = a[...,0]*b[...,0]
    # To comply with the rfft support of multi dimensional arrays
    ar = a.reshape(-1,a.shape[-1])
    br = b.reshape(-1,b.shape[-1])
    cr = c.reshape(-1,c.shape[-1])
    # Note that we cannot use ellipses to achieve that because of 
    # the way `view` work. If there are many dimensions, one should 
    # consider to manually perform the complex multiplication with slices.
    if c.shape[-1] & 0x1: # if odd
        for i in range(len(ar)):
            ac = ar[i,1:].view(np.complex128)
            bc = br[i,1:].view(np.complex128)
            cc = cr[i,1:].view(np.complex128)
            cc[...] = ac*bc
    else:
        for i in range(len(ar)):
            ac = ar[i,1:-1].view(np.complex128)
            bc = br[i,1:-1].view(np.complex128)
            cc = cr[i,1:-1].view(np.complex128)
            cc[...] = ac*bc
        c[...,-1] = a[...,-1]*b[...,-1]
    return c
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

我应该如何将 scipy.fftpack 输出向量相乘? 的相关文章

  • 如何用python脚本控制TP LINK路由器

    我想知道是否有一个工具可以让我连接到路由器并关闭它 然后从 python 脚本重新启动它 我知道如果我写 import os os system ssh l root 192 168 2 1 我可以通过 python 连接到我的路由器 但是
  • 将html数据解析成python列表进行操作

    我正在尝试读取 html 网站并提取其数据 例如 我想查看公司过去 5 年的 EPS 每股收益 基本上 我可以读入它 并且可以使用 BeautifulSoup 或 html2text 创建一个巨大的文本块 然后我想搜索该文件 我一直在使用
  • 使用 Python 从文本中删除非英语单词

    我正在 python 上进行数据清理练习 我正在清理的文本包含我想删除的意大利语单词 我一直在网上搜索是否可以使用像 nltk 这样的工具包在 Python 上执行此操作 例如给出一些文本 Io andiamo to the beach w
  • Pandas 日期时间格式

    是否可以用零后缀表示 pd to datetime 似乎零被删除了 print pd to datetime 2000 07 26 14 21 00 00000 format Y m d H M S f 结果是 2000 07 26 14
  • Python zmq SUB 套接字未接收 MQL5 Zmq PUB 套接字

    我正在尝试在 MQL5 中设置一个 PUB 套接字 并在 Python 中设置一个 SUB 套接字来接收消息 我在 MQL5 中有这个 include
  • 立体太阳图 matplotlib 极坐标图 python

    我正在尝试创建一个与以下类似的简单的立体太阳路径图 http wiki naturalfrequent com wiki Sun Path Diagram http wiki naturalfrequency com wiki Sun Pa
  • 如何在 Python 中解析和比较 ISO 8601 持续时间? [关闭]

    Closed 这个问题正在寻求书籍 工具 软件库等的推荐 不满足堆栈溢出指南 help closed questions 目前不接受答案 我正在寻找一个 Python v2 库 它允许我解析和比较 ISO 8601 持续时间may处于不同单
  • Python 2:SMTPServerDisconnected:连接意外关闭

    我在用 Python 发送电子邮件时遇到一个小问题 me my email address you recipient s email address me email protected cdn cgi l email protectio
  • 从Python中的字典列表中查找特定值

    我的字典列表中有以下数据 data I versicolor 0 Sepal Length 7 9 I setosa 0 I virginica 1 I versicolor 0 I setosa 1 I virginica 0 Sepal
  • 如何使用python在一个文件中写入多行

    如果我知道要写多少行 我就知道如何将多行写入一个文件 但是 当我想写多行时 问题就出现了 但是 我不知道它们会是多少 我正在开发一个应用程序 它从网站上抓取并将结果的链接存储在文本文件中 但是 我们不知道它会回复多少行 我的代码现在如下 r
  • 如何使用 pybrain 黑盒优化训练神经网络来处理监督数据集?

    我玩了一下 pybrain 了解如何生成具有自定义架构的神经网络 并使用反向传播算法将它们训练为监督数据集 然而 我对优化算法以及任务 学习代理和环境的概念感到困惑 例如 我将如何实现一个神经网络 例如 1 以使用 pybrain 遗传算法
  • Cython 和类的构造函数

    我对 Cython 使用默认构造函数有疑问 我的 C 类 Node 如下 Node h class Node public Node std cerr lt lt calling no arg constructor lt lt std e
  • Python3 在 DirectX 游戏中移动鼠标

    我正在尝试构建一个在 DirectX 游戏中执行一些操作的脚本 除了移动鼠标之外 我一切都正常 是否有任何可用的模块可以移动鼠标 适用于 Windows python 3 Thanks I used pynput https pypi or
  • 不同编程语言中的浮点数学

    我知道浮点数学充其量可能是丑陋的 但我想知道是否有人可以解释以下怪癖 在大多数编程语言中 我测试了 0 4 到 0 2 的加法会产生轻微的错误 而 0 4 0 1 0 1 则不会产生错误 两者计算不平等的原因是什么 在各自的编程语言中可以采
  • 仅第一个加载的 Django 站点有效

    我最近向 stackoverflow 提交了一个问题 标题为使用mod wsgi在apache上多次请求后Django无限加载 https stackoverflow com questions 71705909 django infini
  • 根据列 value_counts 过滤数据框(pandas)

    我是第一次尝试熊猫 我有一个包含两列的数据框 user id and string 每个 user id 可能有多个字符串 因此会多次出现在数据帧中 我想从中导出另一个数据框 一个只有那些user ids列出至少有 2 个或更多string
  • Python ImportError:无法导入名称 __init__.py

    我收到此错误 ImportError cannot import name life table from cdc life tables C Users tony OneDrive Documents Retirement retirem
  • 模拟pytest中的异常终止

    我的多线程应用程序遇到了一个错误 主线程的任何异常终止 例如 未捕获的异常或某些信号 都会导致其他线程之一死锁 并阻止进程干净退出 我解决了这个问题 但我想添加一个测试来防止回归 但是 我不知道如何在 pytest 中模拟异常终止 如果我只
  • 更改 Tk 标签小部件中单个单词的颜色

    我想更改 Tkinter 标签小部件中单个单词的字体颜色 我知道可以使用文本小部件来实现与我想要完成的类似的事情 例如使单词 YELLOW 显示为黄色 self text tag config tag yel fg clr yellow s
  • 使用随机放置的 NaN 创建示例 numpy 数组

    出于测试目的 我想创建一个M by Nnumpy 数组与c随机放置的 NaN import numpy as np M 10 N 5 c 15 A np random randn M N A mask np nan 我在创建时遇到问题mas

随机推荐

  • Django的manage.py显示旧命令

    我正在编写自己的 whl 包 在创建了一些新的管理命令并删除了一些旧的命令后 我对自己非常满意 除了在构建我的轮子包之后 与setup py bdist wheel 并将其安装在我的测试服务器上 使用pip install U projec
  • 更改 Drupal 的主题并保留 Garland 作为管理主题?

    如何在不更改管理主题 站点的 contrib 主题和管理界面的 Garland 的情况下将 contrib 主题应用到 Drupal 6 站点 Thanks 转到管理 gt 站点配置 gt 管理主题 在那里 您可以设置管理主题 如果您想对管
  • 使用 inno setup 卸载默认图标

    我正在使用 Inno setup 将卸载图标添加到 开始 菜单文件夹 using the Inno Setup Script Wizard example My program there is a default uninstall ic
  • 无论锁定状态如何写入锁定文件

    有没有办法写入锁定的文件 无论它打开哪个程序 进程 设想 作为服务运行的商业产品会锁定日志文件 服务不能停止 因为这会影响客户 想在文件末尾插入一行作为标记 出现错误 该进程无法访问该文件 因为该文件正在被另一个进程使用 有什么方法可以在锁
  • 如何从 C 文件中获取完整的汇编代码?

    我目前正在尝试找出从相应的 C 源文件生成等效汇编代码的方法 我使用 C 语言已经好几年了 但对汇编语言的经验很少 我能够使用以下命令输出汇编代码 S海湾合作委员会中的选项 然而 生成的汇编代码包含调用指令 这些指令又跳转到另一个函数 例如
  • ActiveRecord find_each 结合 limit 和 order

    我正在尝试使用 ActiveRecord 运行大约 50 000 条记录的查询find each方法 但它似乎忽略了我的其他参数 如下所示 Thing active order created at DESC limit 50000 fin
  • 什么时候在函数外部返回值使用移动与复制?

    读完这篇文章后question https stackoverflow com questions 11914691 copy elision move constructor not called when using ternary e
  • Facebook Open Graph API - og:元标签被忽略

    我正在解决页面标题和图像未包含在 Facebook 点赞中的问题 使用 OG 调试器后 看起来没有og Facebook 正在访问元标签 它说它正在推断og url and og title页面上存在的属性 作为测试 我直接从开发人员文档的
  • 类型错误:无法读取未定义的“想要”属性:

    我一整天都在使用 firebase 成功部署功能 学习如何使用它 我试图看看如果我初始化了另一个部署到同一个项目的目录并且没有任何问题 直到我更新了我的 npm 版本 现在每当我尝试部署时都会收到 发生意外错误 我尝试通过让我自己的用户成为
  • 如何在Listview的OnClickListener中获取被点击项目的位置?

    我想在列表视图上获取选定的列表视图项目 onclick 侦听器 现在我已经实施了onItemClickListener但是当我单击某个项目文本时 它就会被提升 我想在列表视图行单击时引发它 知道如何实现这一目标吗 要获取列表 onClick
  • C# 从 IP 地址字符串中删除尾随“\0”

    我通过 TCP 套接字收到了一个字符串 这个字符串看起来像这样 str 10 100 200 200 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 我应该如何将其解析为 IP 地址 如果我这样做 IPAddress TryPa
  • Spring Cloud Stream 和 RabbitMQ 健康检查

    我有一个使用 Spring Cloud Stream Rabbit 和 Eureka Discovery Client 的简单 Spring Boot 应用程序 该应用程序与 Eureka Server 一起工作正常 并且通过 Rabbit
  • 当 WKExtension.scheduleBackgroundRefresh 应该调用 ScheduledCompletion 处理程序时?

    我正在尝试用这样的行安排后台任务 WKExtension shared scheduleBackgroundRefresh withPreferredDate Date timeIntervalSinceNow TimeInterval 5
  • 我应该向 Proxy.newProxyInstance(...) 提供哪个类加载器?

    我已阅读文档 但我仍然不明白应该提供哪个类加载器作为参数 我尝试了一些选项 但这似乎对编译或代理的行为没有影响 有点令人不安的是 我可以传递任何内容作为类加载器参数 包括null 并且代码仍然可以正常工作 谁能解释一下这一点 并告诉我如果我
  • Android NDK 上“__aeabi_ul2f”的多重定义(libgcc_real.a)

    我收到此错误multiple definition of aeabi ul2f 使用 CMake 在 Android 上编译大型项目时 显然 找到定义的地方之一是 libsmoltcp cpp interface rust a compil
  • 如何将 Spring 应用程序与 Mule ESB 集成

    我想将我的 spring 3 0 应用程序与 Mule ESB Mule3 集成 并为不同的客户端 Net GWT 等 提供这些服务 为了实现这一点 我是否应该将 Spring 应用程序部署为单独的组件并在 Mule 上定义 Endpoin
  • Symfony2.3 在控制器中获取 EntityManager 的更好方法

    我在用着Symfony2 3我目前使用 EntityManager 如下所示 构造 这是使用 EntityManager 的更好方法 构造 或在每个方法中使用 如图所示公共索引Action QuazBar controller class
  • 安装后Eclipse不会显示插件

    我一直在 Ubuntu 9 10 中使用 eclipse 没有任何问题 最近我升级到了10 04 似乎我的 eclipse文件夹被覆盖了 现在我尝试再次安装以前的插件 但是当我重新启动 Eclipse 时 它 们不会出现 我可以在 已安装的
  • 使用多列唯一索引与单散列列

    我有一个表 我需要为多个列提供唯一约束 但是 我还可以基于所有必需字段的散列引入一个额外的列 而不是创建多列唯一索引 那么就数据库性能而言 哪一种会更有效呢 MySQL 建议 https dev mysql com doc refman 5
  • 我应该如何将 scipy.fftpack 输出向量相乘?

    The scipy fftpack rfft函数将 DFT 作为浮点数向量返回 在实数部分和复数部分之间交替 这意味着要与 DFT 相乘 对于卷积 我将必须 手动 进行复杂的乘法 这看起来相当棘手 这一定是人们经常做的事情 我认为 希望有一