稀疏傅里叶变换(sparse FFT)
作者:桂。
时间:2018-01-06 14:00:25
链接:http://www.cnblogs.com/xingshansi/p/8214122.html
前言
对于数字接收来讲,射频域随着带宽的增加,AD、微波、FPGA资源的需求越来越高,但频域开的越宽并不意味着频谱越宽,有限信号内可认为信号在宽开频域稀疏分布,最近较为流行的稀疏FFT(SFFT)是在传统FFT的基础上,利用了信号的稀疏特性,使得计算性能优于FFT。本文简单记录自己的理解。
源代码
一、稀疏FFT
主要是12年MIT的论文:Simple and Practical Algorithm for Sparse Fourier Transform。
核心思想:
SFT 作为这样一种“输出感知”算法,其核心思路是按照一定规则 Γ ( • )将信号频点投入到一组“筐”中(数量为 B,通过滤波器实现 ) . 因频域是稀疏的,各大值点将依很高的概率在各自的“筐”中孤立存在 . 将各“筐”中频点叠加,使 N 点长序列转换为 B点的短序列并作 FFT 运算,根据计算结果,忽略所有不含大值点的“筐”,最后根据对应分“筐”规则,设计重构算法 Γ -1 ( • )恢复出 N 点原始信号频谱。
算法流程:
步骤一:选定一个sigma,实现数据频域重排。
频域重排主要是因为FFT之后,频域大值集中在一起,需要将这一伙打散,保证频域的稀疏性(sparse)。打散之后就能稀疏,依据定理(可见n越大,打散的可能性越大,从这里看,sigma与B还是有关系的,sigma体现了相邻频点的最小间隔,而B决定了每个bucket的宽度):
这里是按一定概率(概率性)将频谱重排,MIT重排思路:
对应代码:
fs = 1000; f0 = 70; N = 256; t = [0:N-1]/fs; x = exp(1j*2*pi*t*f0) ; %%*************Step1:频谱随机重排************** sigma = 19; %inv(sigma) = 27 tao = 3; %permutation perm_num = mod([1:N]*sigma+tao,N)+1; pn = x(perm_num);
当然也有强制(确定性)将数据重排的思路,思路类似,只是实现方式不同:
sigma的选择主要利用性质:
这里sigma^-1是sigma关于模N的逆元(数论倒数)。该点主要说明:变化前后的完备性(信号可重建)。
步骤二:加窗。
这里根据筐的多少(B),平分2pi频域,即带宽2pi/B,理想窗函数为矩形窗,但时域为无限长的sinc函数,需要加窗截断,可选择gausswin。
即win = sinc.*gausswin:
这里不局限于gausswin,满足给定约束的窗函数均可:
步骤三:频域抽取并作FFT
加窗之后,保证了频谱不至于展宽严重,进一步保证了稀疏性。频域降采样:
等价于时域的混叠的傅里叶变换:
故直接在时域进行处理。处理完成后进行FFT运算。
该操作的理论基础为:
不谈加窗操作,x->p->y,可以看出x与y存在映射关系,而y的∑操作可以转化为并行,这样一来可以用多个低速率AD实现一个高速率AD的功能,前提是多个AD需要完全同步。
与上述降速率思想等价的是:如果各AD可做好同步,则多个现有AD的能力,可以做出现有AD难以完成的事。
步骤四:哈希映射
哈希映射的线索为:最终观察的数据Y(k)【即y(n)】->P(k)【即p(n)】->X(k)【即x(n)】:
第一步(y -> p)转化的理论依据:
第二步(p -> x)转化的理论依据:
序列重排后的频域变换为:
这两步解释了算法step4的参数定义:
但实际中并非第一步提到的理想情况,实际情况是:
因此存在一个频谱重建的成功概率问题:
4章节的后半部分主要在证明这个概率问题。
步骤五~六的主要目的,引用原文的说法:
Here are two versions of the inner loop: location loops and estimation loops. Location loops are given a parameter d, and output a set I ⊂ [n] of dkn/B coordinates that contains each large coefficient with“good” probability. Estimation loops are given a set I ⊂ [n] and estimate x I such that each coordinate is estimated well with “good” probability.
步骤五:定位循环
d是一个参数,存在约束
原文取值为:
该步骤的主要操作为:
将 Z ( k )中 dK 个较大值(从大到小依次排列)的坐标( k)归入集合 J 中,通过哈希反映射得到 J 的原像:
这样便得到包含原始频谱坐标的集合,迭代L次:
迭代的目的?
步骤六:估值循环
对于每个k∈Ir,计算频率信息:
步骤5~6主要操作:
至此完成了稀疏FFT(sparse FFT)的整个运算过程,记作sfft_v01。
改进版都是依次大框架进行,不同点主要有三个方面:1)重排的实现思路;2)频谱的重建思路。不再一一展开。
原文示例:
该示例给出了步骤5~6的直观解释:
二、问题记录
仿真验证:n = 256点频信号,sigma = 19,则inv_sigma = 27,假设生成点频信号:
图1为重排的信号,图2为重排加窗(此处窗不够理想,看到长尾,论文中强调加窗,作用在于把尾巴剁掉。),图3为原始信号频谱,图4为降采样之后的频谱。可以看出重排容易引入谐波:以sigma为周期。
原始频谱位置为k=10-1 = 9(下标从0开始),重排后对应频谱位置为mod(sigma*k,n) = 171 = 172-1,与理论分析一致。即原始k点重排后位置:
即时域信号,x(n),重排序号perm_num = mod([1:n]*sigma,n)+1,则对应重排后的傅里叶变换,频谱顺序与原信号频谱X,经过重排序号perm_num1 = mod([1:n]*inv_sigma,n)+1,位置始终相差1(数论倒数的性质决定):mod(sigma*inv_sigma,n)=1
稀疏傅里叶变换(sparse FFT)相关推荐
- 【学习笔记】超简单的快速傅里叶变换(FFT)(含全套证明)
整理的算法模板合集: ACM模板 目录 一.概念概述 二.前置知识 1. 多项式 2. 复数 4. 欧拉公式证明 3. 复数的单位根 / 单位向量 三.FFT 算法概述 四.离散傅里叶变换(DFT) ...
- c++稀疏表sparse table的实现算法(附完整源码)
C++稀疏表sparse table的实现算法 C++稀疏表sparse table的实现算法完整源码(定义,实现,main函数测试) C++稀疏表sparse table的实现算法完整源码(定义,实 ...
- 快速傅里叶变换_计算物理基础:第八章-快速傅里叶变换(FFT)
参考北京师范大学的<计算物理基础> 第八章-快速傅里叶变换 计算物理基础_中国大学MOOC(慕课)www.icourse163.org 1.快速傅里叶变换 1.1 离散傅里叶变换及其变换 ...
- 快速傅里叶变换(FFT)详解
快速傅里叶变换(FFT)详解 (这是我第一次写博,不喜勿喷...) 关于FFT已经听闻已久了,这次终于有机会在Function2的介绍下来了解一下FFT了. 快速傅里叶变换(Fast Fourier ...
- 快速傅里叶变换(FFT)的C#实现及详细注释
快速傅里叶变换(FFT)的C#实现及详细注释 ----------------------------------------------------------------------------- ...
- 11.频域里的卷积——介绍,傅里叶变换和卷积,快速傅里叶变换(FFT)_1
目录 介绍 傅里叶变换和卷积 FFT 介绍 我们将继续讨论频率分析以及如何用频率分量的概念来研究图像.如果你还记得上次我们讲过的基于频率的图像分解的概念.我们通过给你们看这张照片来回忆它(如图).这是 ...
- python的图像傅里叶变换 np.fft.fft2 cv.dft 函数
码字不易,如果对您有所帮助,记着点赞哦! 一. 图像傅里叶变换原理: 原理简介请参考:https://www.cnblogs.com/wojianxin/p/12529809.html 对二维图像进行 ...
- 快速傅里叶变换(FFT)和逆快速傅里叶变换(IFFT)
多项式表示法与卷积 多项式有两种表示方法 系数表示法 点值表示法 系数表示法 就是最普通的表示方法,如 f(x)=a0x0+a1x1+a2x2+......+an−1xn−1f(x) = a_0x^0 ...
- 快速傅里叶变换 (FFT)基础
参考博文 频谱分析是一种将复杂信号分解为较简单信号的技术. 首先来很清楚几个概念: FT(Fourier Transformation): 傅里叶变换.其时域信号.频域信号都是连续的. DTFT(Di ...
最新文章
- DeepLearning:tensorflow 参数初始化和参数保存
- 【译】RAID的概念和RAID对于SQL性能的影响
- java的for循环
- 表单必填标星_怎么用JS做form表单验证,要详细代码,求救!(带星号的是必填项)...
- *【ZOJ - 3703】Happy Programming Contest(带优先级的01背包)
- Emulator: PANIC: Cannot find AVD system path. Please define ANDROID_SDK_ROOT
- php mysql 冒号_php – 使用pdo在搜索变量中使用冒号(:)进行查询
- python项目---数据可视化(02)
- Linux学习笔记-项目部署01
- Phoenix=HBase+SQL,让HBase插上了翅膀
- PLC程序的组成结构
- [转]信息安全相关理论题(四)
- 微信小程序获取外部图标iconfont
- linux下使用PulseAudio获取扬声器的音量和是否静音
- 在matlab中字母的小写转换成大写字母,wps文字怎么将大写字母转换成小写字母
- java输入小写字母_java中怎么实现从对话框输入一个大写字母将其转化为小写字母输出?...
- 苏阳乐队杭州巡演后记
- Vivado 除法器IP核 小数模式(Fractional)下结果的修正
- kafka 中如何保证数据消息不丢失
- 浙江大学计算机柳铮,预告 | 计算机学院第二十次研究生代表大会
热门文章
- Git学习(一)(2015年11月12日)
- 登录表单 参考新浪微博
- vimrc常用配置项
- 每个人都有自己的秘密
- matlab 高分屏 变小,解决Ubuntu高分屏下matlab标题栏(菜单栏)字体过小问题
- 用C语言计算中位数 众数,统计学计算中位数与众数
- php mysql 数据库操作类_php mysql数据库操作类
- java10.0.1怎么安装_关于tomcat:您使用哪个Java? 在Server 2016上安装JDK和JRE(10.0.1); 设置JAVA_HOME和路径...
- C语言编写一个赋值程序,实验2 用C语言编写简单程序——2.1 基本数据处理.doc
- 卷积神经网络(CNN)详解及TensorFlow2代码实现