科学网

 找回密码
  注册
科学网 标签 FFT

tag 标签: FFT

相关帖子

版块 作者 回复/查看 最后发表

没有相关内容

相关日志

[转载]尺度滤波器对的双正交性的检验方法之间的数字误差
SciteJushi 2016-6-12 07:28
原载 http://blog.sina.com.cn/s/blog_729a92140102weau.html 在讨论“复数序列的正交镜像quadraturemirror”等时,曾使用了Matlab-R2011a中的相关函数(序列)的计算程序xcorr。因为,c=xcorr(f1,f2)的结果,直接对应于复数序列f1和f2的平移系之间的内积,所以,尺度滤波器f1与f2的双正交性,是指:在由c的隔2下抽样分得到的两个子序列之中,有一个子序列,只在某单个时间点上有非零实部值(与之相比,虚部和其余时刻上的值,为0或小得可以忽略不计),即实质上为冲激序列。 序列f1和f2与c之间的这种关系,明确了在f1与f2双正交时其频谱函数的一个等式(参见小波名著“十讲”的实系数滤波器组)。据这一频域等式,可以检验和设计FIR或IIR。当容许时域复数序列和非标准直流增益时,还可利用好的频谱函数或FIR,再灵活地搬移f1和f2的频谱。 设W0和W1,是经过了双正交性检验、时间点对齐和增益定标后的已合格的一对尺度滤波器序列。 给它们附加相位频谱(结果对应绝对可和序列,可用有限长噪声样本定义)后,得新滤波器序列 H0、H1,再用正交镜像处理又得其小波滤波器序列分别为G1、G0。然后,对于这4个滤波器,可以再附加另一重相位频谱(《 尺度独立地噪声化相谱的复小波包向量的测试》,2015-06-15),并把结果仍记为H0、H1、G1、G0。它们可构成双正交共轭滤波系统。H0=H1、G0=G1时的特例,是正交共轭滤波系统。居士已在博客中给出了基于Tpwp的大量此类例子,实际上,还反复测验过近百个FIR序列对W0、W1(不等同)的频谱的随机搬移。 图片1.是不用Tpwp而且受限制于实变换时,演示和测验上述方法的一个程序(《给小波工具箱基于FIR的实变换附加多重噪声化相谱》,2016-01-22)。由Matlab-R2011a小波工具箱的函数biorwavf取出其15组双正交尺度滤波器序列,生成了滤波器四件套{H0,H1,G1,G0}。为了直接用于普通实变换的双通道卷积滤波的计算流程,再将H1、G1用其倒序后的结果取代。然后,滤波器对{H0,G0}或{H1,G1}之中,任一对都可作为工具箱内的函数wavedec的输入,按其常规(仅用per模式,注意其时间零点规则)实现离散小波正变换,而另一对,可作为函数waverec的输入以实现逆变换。 图片1.的右上角的图形窗,显示了一组测验结果,证实了高精度可逆变换。其中512条重建误差曲线,分别对应序列长度为2至1024的偶数的被处理信号。最大分解深度可达10。这表明,实变换中的双正交滤波器组的双重噪声化相谱方法正确,也没有使双正交离散小波变换的数字稳定性变得更差。 那个xcorr,是Matlab的信号处理工具箱中的一个函数。近些时候,居士已注意到,数学书中的自相关函数的定义,有时不同于信号处理中的常见概念。在《现代数学手册·随机数学卷》(2000年)第598页中,“相关函数”,似乎对应于信号处理时实际用的“相关系数”,与协方差函数和归一化有关,而且其中二阶矩的表达式中也无序列值的复共轭运算。这将涉及功率谱的估计、傅立叶变换的应用、相关序列与卷积运算的关系等,要求小心均值的问题。不同背景的人们,可能有些不同习惯。又例如,在工科书籍中,不少人用“微分器”、“微分运算”,与“积分器”、“积分运算”对比而言,实际上可能表示数学中的“求微商、求导数函数”。 在Matlab的小波工具箱的文档中,可以见到,将FIR尺度序列倒序后与另一个序列做卷积的方法,例如,“PtheassociatedLagrangeàtrousfilterisasymmetricfilteroflength4N-1.Pisdefinedby P= ...Then,ifwdenotesdbN…,wisasquarerootofP:P=conv(wrev(w),w)wherewisafilteroflength2N.”,使用了函数wrev实现倒序,而正如工具箱计算正交镜像滤波器(函数qmf)一样,它对倒序后的序列都不做复共轭运算。检验复数值的滤波器序列时,要把倒序后的结果取复共轭后,再输入到卷积函数conv。 用FFT和IFFT,可计算两个有限长序列的卷积结果和相关序列,当然就可以直接检验尺度滤波器对的双正交性。 下面,是在Matlab-R2011a中运行的一个例子。从矩阵c中,经过“2800”行的隔行采样的结果,都验证了双正交性。使用conv时的结果(最后一列),可能更精确地刻画以适当整数形式给出的FIR序列。用xcorr或FFT时,在计算过程中引入了无理数和舍入量,可以看到不可避免的数字误差(与最右边的“0”和“2800”比较)的影响。 clear f1= ; f2= ; %f1=f1*exp(3i);f2=f2*exp(7i); c(:,1)=fftshift(ifft(fft(f1,13).*conj(fft(f2,13)))); c(:,2)=xcorr(f1,f2); c(:,3)=conv(f1,conj(f2(end:-1:1))); c c= 2.2864e-0131.6038e-0130 333 -5.2471e-014-1.7109e-0130 -184-184-184 2.1594e-013-3.9252e-0160 158115811581 280028002800 158115811581 1.6884e-013-1.6038e-0130 -184-184-184 -1.9239e-013-5.6288e-0140 333 8.6181e-0143.9252e-0160 新浪赛特居士SciteJushi-2016-06-12。 图片 1.给实变换中双正交滤波器组添加双重噪声化相谱的例子和结果
905 次阅读|0 个评论
[转载]减少DFT算法模式中实变换时的频域乘法的次数
SciteJushi 2016-5-13 08:17
原载 http://blog.sina.com.cn/s/blog_729a92140102wc7i.html 本短文用“实变换时”表明:限于Tpwp中,特别地假定被处理信号和滤波器冲激响应,都是实数序列。不管Tpwp之外的变换,因为使用DFT的条件和方法不一样。这里,不涉及人们从文献中可能见到的奇偶时间点序列分组、FIR以及时间延迟之类的问题,给定的一个被处理序列x和滤波器序列F(高通或低通),必须靠“移位方向参量为空值”另外指明它们都是频域区间 的结果仍等于v。 以博文《减少小波包变换的DFT算法模式中的DFT的次数和点数》(2015-10-26)为基础,减少了部分复数乘法,但是,新增了复共轭、向量倒序的处理内容,改变了那里的“频域数据重复”与频域的乘积运算的顺序。 用这两个函数子程序,覆盖掉以前用的旧子程序,重新运行博文《冗余小波包合成以及降噪试验和能量守恒定标》(2015-12-27)的最后三个降噪实验。可见在Scilab-5.3.3的命令窗中显示了相同的数字(format(‘v’,10) 格式)结果。但是,有些遗憾,没有见到实验的总耗时量明显减少,甚至反而可见高达百分之十几的增加,这可能与引入的其它处理以及系统工作平台中整体的解释、管理等有关。据此类经验,比较算法的性能时,居士希望尽量交代清楚方案、条件、环境。 新浪赛特居士SciteJushi-2016-05-13。 图片 1. 减少了紧凑DFT算法模式中实变换时的频域乘法的次数的程序段 附 . 最近,搜索小波的最新研究进展,在《百度百科》中,注意到了中国已有研究成果“比小波分析创始人Mallat提出的著名Mallat算法快一个数量级以上,而且效果和质量远优于Mallat算法”。看到这很振奋又好奇,然而还未搜到演示资料或方法细述,故特将其复制为留念于下。原来,还有了国际小波分析应用研究中心。虽然,姓名,是居士在博客中多次提到(因个人认为,翻译得已很不错)的小波名著“十讲”的译者,也是前面(2016-04-20)的实变换的滤波器构造中曾用过那篇中文参考文献的作者。但是,居士真不记得当年离开小波的研究时,对这个姓名已有印象,实际上,除去博文已明确提到的部分,至今仍几乎未知其具体研究(注意,在此无褒贬之意)。 : 李建平(电子科技大学教授) 简介?编辑 1986年获重庆大学应用数学理学学士,1989年获西安交通大学软件工程、计算数学双硕士学位,1998年获重庆大学计算机科学理论工学博士学位,1999年至2001年到香港浸会大学做博士后研究,2002年至2003年到法国普罗旺斯大学做访问学者,2005年至2006年到加拿大圭尔夫大学、多伦多大学做访问学者. 教学情况 现为电子科技大学计算机科学与工程学院、示范性软件学院副院长,国际小波分析应用研究中心主任,国际学术进展 International Progress on Wavelet Active Media Technology and Information Processing主编,国际学术期刊International Journal of Wavelet Multiresolution and Information Processing副主编,兼任国家科学技术奖励评审委员、国家自然科学基金项目评审委员、中华人民共和国公安部技术顾问等十几个学术、社会职务,是2004年国际计算机学术大会、第三届国际小波分析及其应用学术大会(ICWAA)、第二届智能体媒介技术国际学术大会(ICAMT)主席,是国际上小波分析与信息处理研究领域十分活跃的知名专家。 科研情况 他在国际上独立提出并系统建立了“小波变换的加速方法”、“矢量积小波变换理论”、“基于小波分析的电子签名系统”等系列理论与方法,该理论比小波分析创始人 Mallat提出的著名Mallat算法快一个数量级以上,而且效果和质量远优于Mallat算法,这为小波理论在信号处理、信息分析等许多方面的应用提供了先进的算法,为小波分析的实时处理和产品化提供了理论基础,为小波分析开辟了广阔的应用前景。他在国际上首次提出“基于‘三大特征’(机器特征、文档特征、人体特征)的信息最安全传输的模型与方法”,为当今研究热点的网络与信息安全做出了重要贡献。他先后主持国家863高技术项目、国家自然科学基金项目等30项。他在国内外著名学术期刊上发表重要论文150篇,被国际三大检索机构SCI、EI、ISTP等检索收录论文48篇,出版学术专著15部,其中2部被多次修订重印。他主持研制的“小波指纹加密系统”、“分布式网络监控系统”等高技术产品产生了广泛的经济效益和社会影响
789 次阅读|1 个评论
[转载]为离散正交变换试凑伪复随机FIR尺度滤波器
SciteJushi 2016-4-20 09:01
原载 http://blog.sina.com.cn/s/blog_729a92140102waam.html 这里称呼“尺度滤波器”,但是不必要求其直流增益为根号2,甚至可能已不该称之为滤波器,而且它是由一系列伪随机数为参数而来,可能正好是通常的高通滤波器,也可能就是单位冲激序列,所以加了“伪”字。 用“试凑”,以声明:靠设想和数字计算的检验,考察复数值FIR序列,未细查有关历史,不要求每步处理的必要意义、约束参数以达到某种性能等等。与这相比,此前研讨中,用的频率响应,更明确些。对于实变换的滤波器构造感兴趣者,可能由《一种新的面向信号处理的小波变换加速算法》(《软件学报》13(07)1338-07,2002)入手。 图片1.的右半部,完整地显示了本文生成FIR滤波器序列(第一个输出变量)的函数子程序TryRandFIR .m,它仅有十几行可执行的普通语句。其6个输入变量I0、SU0、U0、U1、e1、o1依次表示,FIR长度值的二分之一、控制递推“实模”参数(给cos、sin)序列和值的指示、实模参数序列的随机均匀分布的范围、随机“相位”参数序列的三个均匀分布的范围。 图片1.的左半部的命令窗中,显示了用Matlab的xcorr函数计算FIR的自相关序列的快速测验的两组结果。这无需小波变换,方便快速。所有360个FIR序列的检验结果都表明,TryRandFIR的输出序列的偶数平移系,构成了范数为1的正交向量集,当然,这可再与《直面quadrature》(2016-04-06)所述,珠联璧合。 其中第一组测验,表明,当模参数序列和为pi/4、最后3个输入相位参数都为0时,输出序列可能据直流增益为根号2来重定标,合乎了实变换的滤波器构造的旧文献。 第二组测验中,最后的4个参数都为1e6,使随机值的变化范围很大。图片1.的底部的曲线图,显示了记录的8个例子,表明,所实现的FIR序列的幅度谱,已有很大随机性和复杂度。 在本短文的测验中,FIR的长度,可取得2至360的所有偶数。当FIR长度值较大时,其末端的数值,常已极小,实模参数序列不为常数值时的FIR复模值的对数曲线图(semilogy),有拱形轮廓。 用之于离散“小波包变换”的测验程序,PwpRandomFIR.m,如图片2.的右半部所示,用《尺度独立地噪声化相谱的复小波包向量的测试》(2015-06-15)的PwpRandComplexPhas.m改写而成。只修改了附加尺度独立的噪声化相谱之前的那些语句。因为,这里不按尺度滤波器的直流增益定标,直接用FFT、IFFT来生成周期化滤波器对,所以,移去了过去常用的强迫复数化处理、大部分的随机设置。仍保留了设置Wp.RightShift,以随机抽选约60%的变换,用DFT算法模式。在共轭正交镜像(CQM)处理之前,已预添加了一重噪声化相谱。 如图片2.的左半部所示,用了46000多个随机FIR序列的一组结果,证实了高精度离散正交小波包变换。测试涵盖了2至512的所有偶数长度的复随机信号,最大分解深度,已达9。 新浪赛特居士SciteJushi-2016-04-20。 图片 1. 构造复随机FIR滤波器的 程序及其输出序列的自相关测试 图片 2. 用小波包变换测验复随机FIR滤波器的 程序及结果
650 次阅读|0 个评论
[转载]用小波工具箱的实变换试验滤波器组的随机直流分配
SciteJushi 2016-4-12 08:12
原载 http://blog.sina.com.cn/s/blog_729a92140102w9mt.html 被处理信号中的直流成分(均值),完全由尺度(低通)滤波器中输出,然而小波(高通)滤波器,对于直流成分的增益因子为零。即使《用小波包变换测验伪随机频响滤波器组》(2016-03-26),也强令了W0(1)=1,以迫使低通滤波器的幅频响应在零频率处取得最大值。这是小波变换的一种传统印记,也合乎人们对于滤波器的高通、低通之分的习惯,但是,对于某个变换中的完全重建性质和保内性质等,并不是必须的。 用滤波器“组”时,省去“器”字,感到点别扭。用“系统”而省略 “器”时,感觉正常,甚至略好些。称呼中省不省去“共轭”一词,有教材称“正交滤波器组”、“双正交滤波器组”、“低通滤波器H和高通滤波器G组成一滤波器组,使用共轭滤波器组(H*,G*)对原始信号进行分解,然后用(H,G)重构信号”,这些不是致命的麻烦。 一组基向量,在把它的数值都取复共轭后, 结果仍为基,而且在实变换中其实是不变的。给定的一对分解滤波器,也可 被用 为合成滤波器。 离散傅立叶变换、基向量组、小波包变换算法,如水乳交融(当然,这不是看,在程序集中,字符串“fft”多;近期欲了解 “ wavelet”的 进展,可搜索到但是难联接不少网页),这是居士的博客的一个重要标志和开创,已显示了灵活调整滤波器组冲激响应序列的频谱的简便途径。 滤波器组的随机直流分配,指:强行把普通尺度滤波器的直流增益,乘以区间 内的一个随机数,并相应调整了小波滤波器的直流增益,却保持其余频率处的复增益值不受影响。 如果,简单地以尺度滤波器的直流增益为根号2来对整个滤波器组的数据定标,那么,这就须在随机直流分配之前完成。 完整的测验程序mlt_FIR_RandGainDC_DWT.m,以及一组结果,如图片1所示。程序与mlt_FIR_RandPhase_DWT.m(《给小波工具箱基于FIR的实变换附加多重噪声化相谱》,2016-01-22),有基本相同的前半部分。 其中W0是常规正交变换中的FIR尺度滤波器序列,编号1至30,是sym1至sym30;编号31至60,为db3至db32。把W0定标后,变换到频域,并用随机数调整直流增益后,取其IFFT即得F0。把F0倒序,再令奇数指标点的元素数值都反号后,即得F1。 这一对F0、F1,既可用作正变换函数wavedec中的分解滤波器对,也可用为逆变换函数waverec中的合成滤波器对。完全不需要,Matlab工具箱的函数orthfilt、qmf、wfilters或biorfilt来区别生成小波滤波器、分解滤波器、合成滤波器。当输入随机序列x时,可得分解结果C和合成结果x1。 再把F0、F1倒序后,又直接作为waverec和wavedec中的滤波器。以C为输入时的重建结果记为y,以x1为输入时的分解结果记为z。 把x与y的误差(范数比),加上x与z的误差,记入误差矩阵。程序结束时,画出所有误差曲线。 结果表明,使用随机直流分配后的滤波器组时,实正交变换仍具有,完全重建性质,分解与合成是可交换的。 新浪赛特居士SciteJushi-2016-04-12。 图片 1. 用小波工具箱的实变换试验随机直流分配的滤波器组的程序和例子
644 次阅读|0 个评论
[转载]用小波工具箱的实变换测验滤波器组频响的随机组合
SciteJushi 2016-3-16 08:35
原载 http://blog.sina.com.cn/s/blog_729a92140102w6lc.html 已知,两个正交共轭滤波器组A、B的尺度滤波器的幅频响应分别为A(w)、B(w),a、b为两个不同时为零的有限的随机实数。将A(w)、B(w)分别与a、b相乘,并把这二个积的平方和,记为C(w)。称C(w)的平方根,为滤波器组A、B的幅频响应的一个随机组合。 本短文,以Symlets和Daubechies系列滤波器组为A、B的例子,把C(w)作为新的尺度滤波器的幅频响应,附加上多重噪声化相谱,再测验小波变换的重建误差。 完整的试验函数子程序mlt_FIR_RandCombi_DWT.m,以及一组测验了长度为2至512的偶数的随机信号的结果曲线,如图片1.所示。附加噪声化(对称的)相谱和测验重建误差的方法与流程,与《给小波工具箱基于FIR的实变换附加多重噪声化相谱》(2016-01-22)中相同。 第一个已知的尺度滤波器序列,依次为sym1至sym30以及db3至db32。第二个尺度滤波器序列,随被处理信号序列的长度不同,从db1至db31中随机抽取。这两个滤波器序列的长度的最大值,除以2后,即是访问重建误差累加矩阵中的元素(未被访问到就始终为零)的行指针。 幅频响应的随机加权值,在 内均匀分布。 结果证实,由随机组合得到的幅频响应而生成的周期滤波器序列,在工具箱的QMF(orthfilt)处理之前或后再附加了噪声化相谱,它们也仍适用于实正交变换。 如果,使用Tpwp的紧凑DFT算法(《减少小波包变换的DFT算法模式中DFT/IDFT的次数和点数》,2015-10-26),那么,通过在频域乘以单位模序列,可更方便地,随意调整滤波器或被处理信号的相位谱,例如,可取代降噪处理中的循环移位(cycle-spinning)。 经周期化而被延长了滤波器序列,对于使用普通卷积下抽样的经典算法的小波包变换来讲,不是好选择,因为可能带来深度分解与合成计算的耗时量的大增。但是,普通小波变换相当于仅处理到最低频的两个节点的数据(以前曾经用过“即使直接用小波基,分解也需要覆盖过半树”,这类说法是不可取的,尽管,后级分解,须用前级的输出数据的“一半”),所以本文的测验,仍可很快地完成。 新浪赛特居士SciteJushi-2016-03-16。 图片 1. 用小波工具箱的实变换测验滤波器组能谱的随机组合的例子和结果
578 次阅读|0 个评论
[转载]小波无罪Shannon采样定理也没有违法
热度 1 SciteJushi 2015-11-22 08:42
原载 http://blog.sina.com.cn/s/blog_729a92140102vz7w.html 由于某些原因,人们在应用实践中,常未明确考虑基于连续时间小波基函数的被处理数据初始化问题,或无交代,甚至有人误以为Mallat算法就已解决了这一问题。有研究者,对此状况,持较严厉的批评态度,于是有“小波罪恶(wavelet crime)”之说,如图片1的顶部所示。 把离散变换和连续时间信号联系起来,可以更好地理解和开拓应用,但这样常不得不接受近似处理。在与Mallat算法有关的小波变换中,过渡联系的麻烦,表现得尤其突出,难被消除。严格地说,在很多实践中,紧支集连续时间小波基函数,还不算真有用。 习惯于傅立叶分析时,可能想,用矩形窗截取FIR低通滤波器序列的傅立叶变换的一个周期,然后利用频谱的伸缩以获得不同的低通带宽度,就得意味着不同分辨率的子空间簇。然而实际上,在Mallat等的连续时间能量有限信号空间的多分辨分析理论中,尺度函数以及嵌套子空间列,有非直观的特殊的定义方式。人们常难由FIR滤波器组用闭式表达连续时间基函数,只知在极限意义上存在,一个具有某些性质的连续时间尺度函数和小波函数系。 即使,可用解析式表示被处理信号,也不一定容易计算出它与某个紧支集基本函数的平移系的内积,更何况,实际应用中常常只知离散序列 (名著“十讲”第五章后的注释 12,是研究初始化滤波的引子), 难确定其对应的合适函数子空间。这个内积,是连续时间基函数的应用解释的严密性的关键,而也反使数字频率与模拟信号物理频率的精确联系变得麻烦。 如果什么都已知的话,也许 人们 就不需离散小波变换,来做压缩或降噪等事情了。 见到Matlab小波工具箱通过数字计算画个小波函数曲线时,也不清楚其中局部扭结小尖峰,正是好的性质的体现,还是有了较多数字误差。 居士脱离连续时间,直接讨论离散小波包变换,可以使计算理论更独立严谨。这正如,离散傅立叶变换及其快速算法,不必依赖于,连续时间和连续频率概念。由离散时间滤波器组,进入某些函数空间,在其中做理论分析等,是其它的故事。不同的部分,各有价值。 某些问题还未解决好,或未被意识到,或者论述行文失误,都是研发中较常见的事,不意味着“犯罪或邪恶”。这正如,不能说,应用傅立叶变换或Shannon采样定理时的不都定量的近似,是“违法”行为。现代的傅立叶变换理论,已包括很多部分,在数字计算中实际常用的是FFT。各部分连接的关键基础,在长时间的大量文献中 , 同样可能有“不尽人意”之处。当然,批评与自我批评,有事说事,努力解决问题,这是正道。 在傅立叶让他的级数出世后,人们用了相当长的时间,理清傅立叶级数的收敛性问题。 图片1.中关于采样定理的描述,取于国内外的著名教科书或著作,并不是研发初学者之作。 看起来,不算“明确而精辟”。 为什么不用开区间,表示信号的频谱?直观地看:采样使频谱周期延拓时,闭区间有重叠点,使“离散谱”混叠、消失、或不能区分两个边界频率;而且,理想滤波器的频响,在边界频率处不连续,使人感到,不确定、不安全、宜避开。开区间,至少从形式上就排除了麻烦:在边界频率上的正弦信号一个周期内,正好都采到零。 为什么不用半开半闭区间,表示信号的频谱?直观地看:频谱周期延拓时,可无重叠点,也包括了整个频率轴,更像是在频谱的采样离散化时单个周期的两端点中只取一侧。 早期的采样定理中的信号及其傅立叶变换,到底需要什么严格条件? 维基百科,并不排除,定理适用于频率在边界以内的正弦信号(如文后附1.)。“教育部研究生工作办公室推荐研究生教学用书”《数字信号处理》,有专节讨论正弦信号抽样的问题。更一般地讲,定理要面对离散频率成分的问题。 绝对可积函数的傅立叶积分变换,或连续频谱函数,不适用于正弦波。像Dr.Daubechies那样,把带限信号集,视为能量有限 (无穷时间) 信号空间的子空间,陈述采样定理时,也排除了正弦函数等功率有限而能量无限的周期信号。 使用冲激函数串来推演采样定理,这意味着引入了广义函数的傅立叶变换。统一表达了离散谱线和连续区间上的频谱之后,定理就包括了正弦函数。 现在相信,实际信号都是能量有限的,不存在无始无终的永动机式的精确单频率运动过程,人们天天做的采样、数-模转换、滤波都是现实的。然而,信号分析实验中常用的三角函数模型和采样定理中的sinc插值重建函数,都不是物理可实现的。 即使,认为现实中能记录的信号,都既是时限的也是带限的,也不都很精确每个实际问题中的误差,一般也不随Slepian等研究“时-频限”。很多时候,也未必明确所谓的冲激函数串的理想化抽象近似的量度。 Shannon采样理论,用sinc函数做 理想 插值重建(Whittaker–Shannon interpolation formula),违背因果律。很强调了带限,削弱了对现实时限的注意。实际上,采样率变化了,那么,某给定时间区间内的重建信号,受该区间外的样本点影响的程度,也就不同了。假如对插值函数做截断近似,则还需要兼顾函数值的量化误差。 至此,容易注意到,过采样的采样率很高时,可以认为很多“插值函数、重建基函数”都可“更近似于冲激函数”,其差别减小了。这样,使用多分辨分析子空间的小波变换的初始化时,误差也可能减小。 记得当年老师说,凭经验,实用电路中,采样频率一般大于信号成分的最高频率的2.5倍,因为低通模拟滤波器的过渡带不是理想的。田边顽童,不识定理,更无非均匀采样理论,而可能用小石子或谷粒拼出些自以为像人头的图案,指指点点,叫这个狗娃,称那个胖仔。 新浪赛特居士SciteJushi-2015-11-22。 图片 1.小波罪恶之说与采样定理的表述 附 1:Wikipedia中的采样定理(Nyquist–Shannon sampling theorem) 部分的摘录 Sampling is the process of converting a signal (for example, a function of continuous time and/or space) into a numeric sequence (a function of discrete time and/or space). Shannon's version of the theorem states: If a function x(t) contains no frequencies higher than B hertz, it is completely determined by giving its ordinates at a series of points spaced 1/(2 B ) seconds apart. A sufficient sample-rate is therefore 2 B samples/second, or anything larger. Equivalently, for a given sample rate f s , perfect reconstruction is guaranteed possible for a bandlimit B ≤ f s /2. When the bandlimit is too high (or there is no bandlimit), the reconstruction exhibits imperfections known as aliasing. Modern statements of the theorem are sometimes careful to explicitly state that x ( t ) must contain no sinusoidal component at exactly frequency B , or that B must be strictly less than ½ the sample rate. The two thresholds, 2 B and f s /2 are respectively called the Nyquist rate and Nyquist frequency . And respectively, they are attributes of x ( t ) and of the sampling equipment. The condition described by these inequalities is called the Nyquist criterion , or sometimes the Raabe condition . The theorem is also applicable to functions of other domains, such as space, in the case of a digitized image. The only change, in the case of other domains, is the units of measure applied to t , f s , and B .
1129 次阅读|2 个评论
[转载]减少小波包变换的DFT算法模式中DFT/IDFT的次数和点数
SciteJushi 2015-10-26 09:33
原载 http://blog.sina.com.cn/s/blog_729a92140102vwr9.html 在小波(包)和滤波器的普通理论中,常见的是Z变换和傅立叶积分变换。为应用方便,基于子带分解塔、周期化共轭滤波器组和Wickerhauser等的最佳基思想,居士建立有限长复向量空间的小波包分割与覆盖以及多分辨分析。其重要意义,不受限制于,连续时间小波函数和无穷时间序列滤波。在博客中,已包括了许多创新,不明现连续时间的基函数(多分辨分析和小波的纯粹数学理论的主要构造目标)、Z变换域分析以及滤波器阵的反混叠条件(anti-aliasing condition)等内容,而反复强调了关键词“基向量”。因为其尺度、时间和频率都是有限的离散量,所以这种小波包变换,与离散傅立叶变换的联系,更直接,工整自然,当然,其理论基础需要Z变换和连续频率。 以《随机设置小波包变换的DFT算法模式》(2015-05-13)为基础,还可使离散傅立叶变换(从通常的数字信号处理教科书即可知其所需性质),与小波包变换的算法实现之间,更紧密地结合。例如: (1).底层(相应于常规的“卷积、滤波、下抽样、插0”处理)的子程序可输入和输出频域数据,供调用程序缓存和重复利用,以免相邻尺度的运算过程中的反复的DFT和IDFT处理。尤其在信号重建中,一般并不需要中间经过的节点上的时域数据。高通支路输出与低通支路输出的相加,在时域或频域均可实现。 (2).用频域半周期折叠平均,代替时域隔2下抽样 ( 2- 抽样,二个点中取一舍一) ,使获得时域系数的IDFT的点数减半,回避下抽样时要抛弃的一半。 (3).用频域数据“重复”,代替时域插0变抽样率,绕过了插零,可使求频域系数的DFT的点数减半。 新旧两套方案对照,如图片1.的右半部的同名文件所示。 这种调整,不可避免地要改变数字误差,就不妨顺便用矩阵乘法代替用求和函数sum做非DFT模式的分解。调用程序,需较大改动,变得更复杂些,尤其在计算所谓的双向树时,才能兼容地管理时域数据和频域数据以及非DFT模式,但是,这不影响此前讨论的应用和测试程序运行及判断。 重新运行《用小波包变换试验多孔滤波器组的零频点移位》(2015-09-09)和《用双重噪声化相谱的多孔滤波器组测验小波包双向树》(2015-09-25)的测验程序,其状况如图片1.左上半部所示。 第二测验程序的耗时仍很长,其中约40%的测试不使用DFT算法模式。仅改动其中的语句“if abs(Wp.RightShift)5e5, Wp.RightShift= ;”即强迫全部用DFT算法模式处理,可重命名函数和运行测试,则耗时从2286秒降至了568秒。 新浪赛特居士SciteJushi-2015-10-26。 图片 1. 更紧凑DFT算法模式及其三个测验程序的情况
736 次阅读|0 个评论
[转载]从零化相谱的Shannon小波包分解树上取FFT复谱数据
SciteJushi 2015-7-15 08:19
原载 http://blog.sina.com.cn/s/blog_729a92140102vqn4.html 可以使小波包变换的滤波器序列的相谱噪声化,这表明, 相位谱可能有多种意义特殊的设计选择 。例如,可能在使用Shannon复小波包时,干脆把相位谱全变为零序列,即:按Tpwp的DFT算法模式中的滤波器组的准备,把普通小波包向量组用FFT变换到频域;在频域只保留复模序列,而抛弃所有非零相位;需要用移位基向量组做内积的分解算法时,就再由IFFT把模序列变换回时域。 这种“零化相谱”处理,把Shannon小波包的基本序列的频谱(周期化滤波器组的频率响应)转化为非负的实数序列。在做小波包分解时,最深分解尺度上的小波包系数,乘以信号序列长度的平方根后,就可能等于用工作平台的FFT函数求得的信号的复数谱数据。 离散傅立叶变换,对人类的幸福和文明的进步,做出了极重要的贡献。居士不愿想象,如果没有离散傅立叶变换,那么科技发展现在处于何种落后状况。小波包变换,与离散傅立叶变换,和谐结节,这令人愉快,似乎平缓了草丛里的长蛇、夜幕中的蚊鸣的惊扰。 风雨过,云天蓝;瓜果熟,稻穗出;人烟稀,宅院空;山林间,新坟旧冢,点点安详。 在小波的理论和应用的通常文献中,正交变换存在著名的非线性相位的问题,由FIR低通滤波器转换到高通滤波器时还有时移,由此,人们就可能怀疑:把周期化的低通滤波器和高通滤波器序列的相位谱,都置为零,而且允许零频率点偏移向低频通带的边缘,如此反常又极端的措施,是可行的。观念或程序软件,局限于FIR滤波器组,且用其支集长度限定了分解的深度,也就难见到在小波包系数集之中的信号频谱序列。 测验程序PwpComplexP0ShanSpec.m,由《从噪声化相谱的Shannon小波包分解树上取FFT幅谱数据》(2015-07-11)中的PwpComplexShanSpec.m改写而来。如图片1.所示,有两处改动,一是做上述的相谱全零化的处理;二是,在比较最深尺度的小波包系数与FFT结果的误差时,不再需“取模”的运算,然而直接比较“复数序列”间的误差。 虽然不能定义虚部不为零的复数的大小,但是,用Matlab内置的排序函数子程序sort,可以处理复数序列。模不同的复数的顺序,由模确定;在模相同时,才利用相角,进一步确定顺序。 图片1.中的Matlab命令窗内,显示了用大量的复随机信号(程序输入200、88)的测试的结果。重建误差、复谱序列的误差,仍都极小,小于万亿分之一。 新浪赛特居士SciteJushi-2015-07-15。 图片 1. 零化相谱的Shannon小波包分解树上的FFT复谱数据的测验程序和结果
793 次阅读|0 个评论
[转载]从噪声化相谱的Shannon小波包分解树上取FFT幅谱数据
SciteJushi 2015-7-11 17:08
原载 http://blog.sina.com.cn/s/blog_729a92140102vqjs.html 使用“取FFT幅谱数据”一语,呼应博文《用基向量从随机设置小波包变换树上捕捉DFT负频率的幽灵》(2014-03-27)。给出特殊例子,最深分解尺度上的小波包系数,乘以信号序列长度的平方根,则其模序列可等于用工作平台的FFT函数求得的幅度谱数据。 只涉及离散Shannon复小波包(《试造离散Shannon复小波包变换》,2015-07-05),但是,仍允许通带中心随机偏移,而且可用《尺度独立地噪声化相谱的复小波包向量的测试》(2015-06-15)中的方法,使小波包向量组的相位谱噪声化。仅做正交变换,把滤波器组的直流增益的模值,固定为根号2。 测验程序PwpComplexShanSpec.m,如图片1.所示。有两个输入参数,依次为:指定一个信号长度后,将进行随机测试的次数;随机数种子。有三个输出矩阵,依次为:用最深尺度的小波包系数所重建的信号,与被分解信号之间的误差;将被分解信号,做FFT,然后直接用IFFT所重建的信号的误差;最深尺度的小波包系数的复模序列,乘以信号长度值的平方根后,与直接由FFT输出所得的幅度谱序列之间的误差。输出矩阵的行号,即是随机测试的序号;其列号,为可行的最大分解深度,对应了信号长度。 特地插入“错误”报告语句error(‘Orthogonal ?’),强调“是正交变换”。 计算幅度谱数据的误差之前,将两个模值序列做排序处理。 预定,约40%的测验用DFT算法模式,用左移位向量组或 右 移位向量组做内积的分解方法的信号,各约占30%;被处理信号长度为2的整数次幂;幂指数,即是最大分解深度,取为1至12。 图片1.的右上部的命令窗中,显示了大量随机测试(程序输入200、88)的结果,所有误差都极小(小于万亿分之一,1e-12)。 新浪赛特居士SciteJushi-2015-07-11。 图片 1. 噪声化相谱的Shannon小波包分解树上的FFT幅谱数据的测验程序和结果
815 次阅读|0 个评论
[转载]随机设置小波包变换的DFT算法模式
SciteJushi 2015-5-13 08:02
原载 http://blog.sina.com.cn/s/blog_729a92140102vngc.html 此前,在讨论小包分解与合成计算时,用变量Wp.RightShift控制生成基向量组时移位的左、右方向两种选择。在Tpwp的函数子程序( 如图片1.的右上角所示 )中,添加几个语句 (即if isempty(Right_Shift)…end部分),就直接增设了DFT(离散傅立叶变换)算法模式。用工作平台自带的计算DFT和逆变换IDFT的函数fft和ifft。因为,不必究等效地实现了以前的某部分计算,也不强调其对应外部的反转、移位关系等,所以居士列其为小波包变换的另一种设置选项。 虽然,使用的仍都是已知的时域尺度序列数据,但是,Tpwp明显不排除:直接解算和构造滤波器组的“离散频率”序列。 反复用“随机设置”一词,针对深奥呆板的形象。在不同的编程语言环境中,也不必精确一一对等地,实现那些所谓的设置。不证明某个选择特别正当,但是,始终以Matlab-R2011a工具箱中“实数域的普通离散小波变换”的“标准模式和算法”的重建误差为参照,表明它们都可行。某个选择可能只构成另一些设置的等效计算,也可能意味着用了不同的小波包基,可能使“小波包林”扩张了。 可用快速傅立(或付、里)叶变换FFT计算卷积,这在信号处理教科书中是常识。在小波变换理论与应用文献中,离散信号的分解与重构,都使用了线性卷积滤波器。自然地,很早就有研究论文,试用FFT来计算离散小波变换。本短文更表明, Tpwp的演绎机制的特点,使增设DFT算法模式,尤其简捷。 小波工具箱中的“数据边界延拓模式”以及Wavelab850中的“延拓后用函数filter”的方法,把一个圆卷积问题,也反而转化成了线性卷积来处理。然而,Tpwp理论把线性卷积问题首先周期化,再推演有限维酉空间的小波包分割与覆盖,所以其计算与FFT更靠近。 对被分解序列做FFT,再与基向量的频域序列F相乘,把乘积序列的逆傅立叶变换IFFT后(注意,也已算出要抛弃的部分)的下抽象样值,作为分解结果。对输入序列补插零值后,做FFT,再与基向量的频域序列相乘,把乘积序列的IFFT输出,作为合成结果。 验证方法的程序之一, PwpRandComplexFFT.m ,由《随机设置小波包变换的强迫复数化的测试》( 2014-11-25)中的 PwpRandComplex.m 改写而成。从图片1.的右下角,可见主要改变,在于插入了条件语句 if abs(Wp.RightShift)5e5, % about 60% tests Wp.RightShift= =TpwpGHfft(... phcoef,pgcoef,rphcoef,rpgcoef,Wp.DatL,Wp.DecL); end % FFT, on all the basis vectors 其中,把 Wp.RightShift,强行改为空矩阵作为标志,指示按 DFT算法模式使用基向量组,否者,仍用旧移位方式;子程序TpwpGHfft用FFT,预先将基本基向量组变换到频域。这是正确使用DFT算法模式,需要做的两件事情。其余程序语句,与 PwpRandComplex.m比, 不变。 一般说来,DFT要求整个处理系统用复数域,而且,函数子程序集,要能兼容控制复数化共轭滤波器组。若在降噪以及预置基的函数子程序中增加TpwpGHfft控制,那么,在降噪和分类识别等试验中,就可使用DFT算法模式了。新函数子程序包TpwpSubs.bin,比旧包大了11KB;TpwpSubs_E01.bin,约增大了1KB。仍用完全安装的Scilab-5.3.3制作,不从互联网下载任何工具箱和优化库文件,也无需用户库。 重编译 子程序集时,尽量少改,即使对于曾测到 rigrsure应用中的被零除的状况也保留(系数加权时可滤掉)。 用高级语言 编 程,且一般做“概念性的、演示性的”实现,不必需移位、数据寻址等的效率优化。 但是,顺便,遵循 Matlab编辑器的建议,在某些子程序中,把原未经初始化然而在循环语句中却需不断改变大小的向量,首先初始化了,这可能影响旧试验的 总耗时量。 生成图片1.中的曲线图,总的耗时约为1955秒,验证了兼容性、双向重建性质、双向保内性质。使用了128组尺度序列、2560组随机复数化滤波器、2560次深度为8的完整树小波包分解、5120次约束树分解、10240次小波包合成。每次分解的复随机信号以及用于合成的复随机系数序列的长度,随机地取为1、3、5、7、9或11与256的乘积。其中约百分之六十,用到了DFT算法模式。FFTW库的设置为estimate,这是Matlab默认的(其帮助文档说“The default method is estimate…If you set the method to 'estimate', the FFTW library does not use run-time tuning to select the algorithms. The resulting algorithms might not be optimal.”)。 用同样的输入参数,运行原 PwpRandComplex.m ,即无 DFT算法模式时,总的耗时约为:4381秒。 在《续用小规模人工神经网络试做个对留数不敏感的分类器》(2015-04-02)所用的测验程序WpdcNN1.m中,把参数Wp.RightShift=1改为Wp.RightShift= ,即将“右移位方式”改为“DFT算法模式”后,再运行,可见降噪性能非常接近然而可节约约百分之二十五的时间。 被处理信号和滤波器组都是实数值序列时,如果不引入DFT,那么分解与合成的结果,也都自动为实数值。然而,经复变换过程后,数字误差可能引入复数,尽管虚部可以很小。虚部值依赖于机器系统精度,也与实现傅里叶变换的子程序对对称性、虚部、实数性等的识别和控制有关(据居士已有经验,Matlab做得很好),这可能使兼容性问题和应用复杂化,总的处理时间也可能反而增加。这时,就不宜选DFT算法模式 。 涉及fft和ifft时,不用Scilab-5.5.0。若用,就除去其FFTW模块的库文件libfftw3-3.dll (资源管理器显示,2311KB, 2013-8-30)后,再重启动,使用默认的标准FFT (standard fft by default),才可实现类似于图片1.中那样的测验。这个文件,约为在5.3.3版中旧库(805KB, 2009-10-29)的3倍大,可能有了相当充分的优化控制,而与居士所用的机器系统难兼容。 新浪赛特居士SciteJushi-2015-05-13。 图片 1.DFT算法模式以及用强迫复数化滤波器组的测证结果(Matlab版)
985 次阅读|0 个评论
学三角多项式逼近和快速傅立叶变换(一)
热度 2 cambaluc 2015-4-8 17:02
《数值分析》中有这么一节内容,简单了解一下基本原理,会用matlab调用fft(),多数同学也能学会。但如何学好学清楚这节内容呢,我看如下几条很关键。 (1、要认识到算法的重要性,应用的广泛性,这个不用多讲,学的时候一定要认识到学好它的必要性,不然学的时候模糊,将来用的时候就不顺手。 (2、得亲自找几个数手工算一算,亲自编写一下算法,这是最好的学习方法,也是验证是否真正掌握算法、编写的算法是否有效的方法。 引用一句话——“人们常常说,除非一个人能将某事教给其他人,他才算是真正地了解了这件事。实际上,直到一个人可以将某事教授给计算机,他才算真正地了解了这件事。---Donald Knuth 1974”;所以,照着算法自己编一下程序,是最好的学习方法。 (3、把学习这个算法当作一件有趣好玩的事来做。曾看过一段国外视频公开课,那个老外说:“它真的如此美妙,所以我说这好比一个人生命中,并不一定要有傅里叶算法,但我非常欣赏它,以及它背后的所有思想。” 以下我摘录整理一下我的讲稿和程序,希望对初学者有所借鉴。 ==================================================== 1. 首先要理解傅立叶级数、傅立叶变换的数学表达式。多看相关数学书,虽然不同的书不同作者可能用不同的符号和形式表示,初学时容易引起混淆,但多看、多做比较,可以更清楚地掌握数学语言所表达的意义。(多到图书馆翻翻书) 2. 对傅立叶变换的应用要搞清楚,学任何算法,知道问题的本质,就算解决了一半。傅立叶变换在信号处理中实现时间域=频率域 的转换,是实现滤波的基础算法;在函数逼近中,可把函数展开成傅里叶级数,实现用简单的三角函数替代复杂或难表示的函数,也是进行近似计算的基础算法。傅立叶变换的应用很广泛,网上也有很多直观有趣的图示(如下面两张图)。 3. 在此阐述一下,函数f(x)的傅立叶级数系数的表示和计算。 我看可以从三个方面来理解,这样便于记忆。(当然先要知道正交函数的妙处,正交函数内积为0可简化很多事。) (1) 高数教材中一般都先入为主地告诉你,函数可展成三角级数,系数an怎么定呢?用cosnx乘两端,再从-π到π逐项积分,根据三角函数系的正交性,其余各项为0,得系数an。单从这外角度来看,多数同学虽然可理解,但总感觉太“数学”了,学了《数值分析》还可从以下函数逼近角度理解。 (2) 最佳平方逼近 用正交函数族作最佳平方逼近:求系数a组成的多元函数极值,得法方程,因所选的基函数为正交函数,所以解可直接写为: ak=( f(x), φk(x) ) / ( φk(x), φk(x) ), k=0,1…n. 根据内积的定义直接写了傅立叶系数的表达式。这样更好地理解傅立叶级数。 (3) 函数向量投影到基向量 有时从几何的角度来理解这一问题更直观。把函数作为一个向量,投影到基函数(基向量)上,傅立叶系数就是在各基函数向量上的长度坐标(射影),也就是坐标值。形式如:prjba=|a|cos∠(a,b)=a.b/|b|。函数向量f(x),投射到函数向量cos(nx)上,an=(f(x),cos(nx))/(cos(nx),cos(nx)),够直观吧。(? 4. 离散傅立叶变换(DFT) 学习快速傅立叶变换,当然先要明白普通的离散傅立叶变换程序的编写。 (1) 复数域表示 复数的有限离散形式 (2) C语言程序 先定义数据结构 typedef struct{ double R; double I; }COMPLEX; COMPLEX add(COMPLEX a,COMPLEX b) { COMPLEX c; c.R=a.R+b.R; c.I=a.I+b.I; return c; } COMPLEX sub(COMPLEX a,COMPLEX b) { COMPLEX c; c.R=a.R-b.R; c.I=a.I-b.I; return c; } COMPLEX mul(COMPLEX a,COMPLEX b) { COMPLEX c; c.R=a.R*b.R-a.I*b.I; c.I=a.R*b.I+b.R*a.I; return c; } // 定义变量及数组,分配内存,参数传递见下文程序 (文件太大,转下一篇)
个人分类: 算法|13945 次阅读|3 个评论
地震波的频谱(FFT)
热度 2 Zhouwenju 2013-9-28 20:55
从fourier transform 开始: (1)time to frequency (2)frequency to time 把(2)用复数表达: S(f) = a(f) + ib(f) 求S(f) 的模(a^2+b^2)^1/2 得到振幅谱 求S(f) 的幅角 arctan(b/a) 得到相位谱 实现: SAC FFT命令 unwrap 命令
个人分类: 基础知识|1060 次阅读|3 个评论
遗忘了的,FFT的一个小角落
moset 2012-11-13 00:22
遗忘了的,FFT的一个小角落
快速 傅立叶变换( FFT) 中最常见的基 2 算法,有着蝶形的结构。基 2 的变换序数有着倒序的特点。所说的倒序,就是序数的 2 进制表达,变换后的序数恰好是其表达的中心镜像对称。如 16 点的基 2 变换 ,序数“ 11 ”二进制 1011 ,那么变换后的序数在 11 这个位置就是“ 13 ”二进制 1101 。 蝶状图的序数的算法方式:一级分为奇数偶数两组,偶数组在上方,奇数组在下方。二级把原先分好的偶数组除 2 ,再分奇数组与偶数组,原先的奇数组减一后除 2 ,再分奇数组,偶数组。依次类推,直到分组无法再分。应用 FFT 许多年了,却没认真考虑过为什么它的序数最后是二进制的倒序。 以 16 点 2 基为例 一级: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 二级: 0 2 4 6 8 10 12 14 1 3 5 7 9 11 13 15 三级: 0 4 8 12 2 6 10 14 1 5 9 13 3 7 11 15 四级: 0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15 一级与四级二进制表达分别: 0 : 0000 0 : 0000 1 : 0001 8 : 1000 2 : 0010 4 : 0100 3 : 0011 12 : 1100 4 : 0100 2 : 0010 5 : 0101 10 : 1010 6 : 0110 6 : 0110 7 : 0111 14 : 1110 8 : 1000 1 : 0001 9 : 1001 9 : 1001 10 : 1010 5 : 0101 11 : 1011 13 : 1101 12 : 1100 3 : 0011 13 : 1101 11 : 1011 14 : 1110 7 : 0111 15 : 1111 15 : 1111 上述的序数的算法结构,用 MATLAB 语言描述大约如下: function h=ap(h,n) % n: 2^n for i=0:n-1 T=2^(n-i); for j=0:(2^i-1) h((1+j*T):(j*T+2^(n-i)))=ac(h,1+j*T,n-i); end end end function tmp=ac(h,start,n) tmp=zeros(1,2^n); for i=1:2^(n-1) tmp(i)=h(start+2*(i-1)); tmp(i+2^(n-1))=h(start+2*i-1); end end h=0:15; ap(h,4) ans = 0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15 倒序方式的MATLAB程序: N=4; p(1:2^N,N:-1:1)=dec2bin(0:(2^N-1)); re=bin2dec(char(p)); re = 0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15 昨儿,莫名想到这个倒序序数的事(也许是 11.11 二进制形式和蝴蝶的各种联想吧)。之后动笔排了好些次,也就得到为什么是倒序的原因了。 C++ 程序描述倒序算法如下: int sequ( int *re, unsigned int n ) { unsigned int i; int x = 0 ; do { for ( i = 0;i n; i++ ) { re |= ( ( x i ) 0x01 ) ( n - i - 1 ); } } while ( ++x ( 1 n ) ); return 0; } 当一个序数 x ,按奇偶分组( x16 ): 最后一级,相对于频域的序数的位置是红还是黑,取决于 (x3)1 ,也就是变换后位置 (H 0)1 ,第三级取决于 (x2)1 ,也就是变换后位置 (H 1)1 ;第二级取决于 (x1)1 ,也就是变换后的位置 (H 2 )1 ;第一级是红是黑取决于 (x0)1 ,也就是变换后的位置 (H 3)1 。同样的,其它 2 基变换的序数也满足这样的规律,由 0 到 2^n-1 的顺序序数依次奇偶分组后的序数必然是倒序。 以此文悼念那些业已失去的,浑浑噩噩的日子。
4799 次阅读|0 个评论
MIT研究出比FFT更快的傅立叶变换算法
热度 7 毛宁波 2012-1-19 08:44
MIT研究出比FFT更快的傅立叶变换算法
据MIT新闻网站报道,MIT的科学家研究出比FFT更快的傅立叶变换算法。 The Fourier transform is one of the most fundamental concepts in the information sciences. It’s a method for representing an irregular signal — such as the voltage fluctuations in the wire that connects an MP3 player to a loudspeaker — as a combination of pure frequencies. It’s universal in signal processing, but it can also be used to compress image and audio files, solve differential equations and price stock options, among other things. The reason the Fourier transform is so prevalent is an algorithm called the fast Fourier transform (FFT), devised in the mid-1960s, which made it practical to calculate Fourier transforms on the fly. Ever since the FFT was proposed, however, people have wondered whether an even faster algorithm could be found. At the Association for Computing Machinery’s Symposium on Discrete Algorithms (SODA) this week, a group of MIT researchers will present a new algorithm that, in a large range of practically important cases, improves on the fast Fourier transform. Under some circumstances, the improvement can be dramatic — a tenfold increase in speed. The new algorithm could be particularly useful for image compression, enabling, say, smartphones to wirelessly transmit large video files without draining their batteries or consuming their monthly bandwidth allotments. Like the FFT, the new algorithm works on digital signals. A digital signal is just a series of numbers — discrete samples of an analog signal, such as the sound of a musical instrument. The FFT takes a digital signal containing a certain number of samples and expresses it as the weighted sum of an equivalent number of frequencies. “Weighted” means that some of those frequencies count more toward the total than others. Indeed, many of the frequencies may have such low weights that they can be safely disregarded. That’s why the Fourier transform is useful for compression. An eight-by-eight block of pixels can be thought of as a 64-sample signal, and thus as the sum of 64 different frequencies. But as the researchers point out in their new paper, empirical studies show that on average, 57 of those frequencies can be discarded with minimal loss of image quality. Heavyweight division Signals whose Fourier transforms include a relatively small number of heavily weighted frequencies are called “sparse.” The new algorithm determines the weights of a signal’s most heavily weighted frequencies; the sparser the signal, the greater the speedup the algorithm provides. Indeed, if the signal is sparse enough, the algorithm can simply sample it randomly rather than reading it in its entirety. “In nature, most of the normal signals are sparse,” says Dina Katabi, one of the developers of the new algorithm. Consider, for instance, a recording of a piece of chamber music: The composite signal consists of only a few instruments each playing only one note at a time. A recording, on the other hand, of all possible instruments each playing all possible notes at once wouldn’t be sparse — but neither would it be a signal that anyone cares about. The new algorithm — which associate professor Katabi and professor Piotr Indyk, both of MIT’s Computer Science and Artificial Intelligence Laboratory (CSAIL), developed together with their students Eric Price and Haitham Hassanieh — relies on two key ideas. The first is to divide a signal into narrower slices of bandwidth, sized so that a slice will generally contain only one frequency with a heavy weight. In signal processing, the basic tool for isolating particular frequencies is a filter. But filters tend to have blurry boundaries: One range of frequencies will pass through the filter more or less intact; frequencies just outside that range will be somewhat attenuated; frequencies outside that range will be attenuated still more; and so on, until you reach the frequencies that are filtered out almost perfectly. If it so happens that the one frequency with a heavy weight is at the edge of the filter, however, it could end up so attenuated that it can’t be identified. So the researchers’ first contribution was to find a computationally efficient way to combine filters so that they overlap, ensuring that no frequencies inside the target range will be unduly attenuated, but that the boundaries between slices of spectrum are still fairly sharp. Zeroing in Once they’ve isolated a slice of spectrum, however, the researchers still have to identify the most heavily weighted frequency in that slice. In the SODA paper, they do this by repeatedly cutting the slice of spectrum into smaller pieces and keeping only those in which most of the signal power is concentrated. But in an as-yet-unpublished paper , they describe a much more efficient technique, which borrows a signal-processing strategy from 4G cellular networks. Frequencies are generally represented as up-and-down squiggles, but they can also be though of as oscillations; by sampling the same slice of bandwidth at different times, the researchers can determine where the dominant frequency is in its oscillatory cycle. Two University of Michigan researchers — Anna Gilbert, a professor of mathematics, and Martin Strauss, an associate professor of mathematics and of electrical engineering and computer science — had previously proposed an algorithm that improved on the FFT for very sparse signals. “Some of the previous work, including my own with Anna Gilbert and so on, would improve upon the fast Fourier transform algorithm, but only if the sparsity k” — the number of heavily weighted frequencies — “was considerably smaller than the input size n,” Strauss says. The MIT researchers’ algorithm, however, “greatly expands the number of circumstances where one can beat the traditional FFT,” Strauss says. “Even if that number k is starting to get close to n — to all of them being important — this algorithm still gives some improvement over FFT.” 引自: http://web.mit.edu/newsoffice/2012/faster-fourier-transforms-0118.html
个人分类: 其他|7386 次阅读|8 个评论
ubuntu下快速傅立叶变换(FFT)的fortran库的建立
zhoufcumt 2011-11-8 10:47
因为要对GPS时间序列做频谱分析,用到FFT,因此试图建立FFT的库来调用。 下载双精度的fft源代码 www.netlib.org/fftpack/dp.tgz 解压后,进入解压目录,修改Makefile: FC=g77 改为 FC=gfortran install: lib$(LIB).a mv lib$(LIB).a /usr/local/lib rm *.o 去掉rm *.o installshared:lib$(LIB).so mv lib$(LIB).so /usr/local/lib rm *.o ldconfig 去掉rm *.o 改完后保存,make 生成 libdfftpack.a, 接着 make shared 生成libdfftpack.so, make install 将libdfftpack.a拷贝到/usr/local/lib下,make installshared 将libdfftpack.so拷贝到/usr/local/lib下,make test测试安装是否成功, make clean即ok。 具体用法现在还没熟悉,感觉应该和lapack库类似,按照他封装的函数格式调用即可,可参照doc了解下函数结构,我会尽快将体验情况更新到博客里。
个人分类: Fortran|8019 次阅读|0 个评论
快速Fourier变换(FFT)与BEM加速
热度 2 xiaojy 2010-10-12 13:16
用快速Fourier变换(FFT)来加速BEM,这应该是我见到的最简单,但又不失高效的一种快速边界元方法了。不需要复杂的数据结构,与其他快速方法具有差不多的计算效率(当然与处理的问题有一定关系)。 我了解一些,希望与有兴趣的同行讨论。
个人分类: Publications|4461 次阅读|5 个评论
复信号的谱分析
热度 4 huarong1940 2010-3-16 00:01
谱分析技术 在光学、声学、电学、力学和信息技术(声像处理、数据压缩)等科技领域里发挥着重要的作用,这是众所周知的。通常的的谱分析技术,大多建立在 实函数 的 Fourier (傅里叶)变换的数学理论基础之上,即使图像处理用到二重 Fourier 变换,也仅以实函数为输入。然而,经典的 Fourier 变换式是 复指数展开 公式,原本是描述复函数及其双边频谱之间的转换, 变换对的双方皆为复数 。常见的谱分析对象为一维的实信号,一维(实)函数是二维(复)函数的特例,它经 Fourier 变换所得的双边频谱是对称的。因此,古老的谱分析技术只利用了 Fourier 变换的单边频谱,仅使用 Fourier 的 三角级数展开 式即可。人们熟知:单边谱的每一频率成分是 谐波函数 ;但几乎不知道: 双边谱的每一频率成分是旋转矢量 。(参阅文后的上传插图《Fourier变换复指数展开式的几何意义》) 上世纪八十年代初,我在对转轴的回转精度的研究中,认识到转子的径向误差是一个随时间变化的二维的复向量,其检测结果是一个二维的复信号,将该复信号的数据直接输入 FFT ,所获得的输出是一个非对称的双边频谱,同时弄明白了双边谱的几何意义。由此,开创了 复信号谱分析 的研究方向,并取得一系列的研究成果,如 回转精度理论 、 机构综合的谱分析方法 、 极形轨迹发生器 等。 Fourier 变换原本是《信号分析》学科的重要理论基础,但是,许多从事《信号分析》研究和教学的学者,至今并未全面、深入地理解堪称经典的 Fourier 变换。迄今,他们仍认定“ 复信号是实际不存在的 ”、 双边频谱的负半边“既无几何意义,又无物理意义 ” 。甚至还把这陈旧的错误论点写进教科书。其错误的根源在于,将研究分析的信号局限于一维空间之内。 许多理工科学者自以为掌握了处理一维、二维甚至多维空间的数学理论,但是一联系到实际,却只习惯于一维的思维和计算。 为了全面理解 Fourier 变换,仅须把思维扩展到二维空间。 扩展思维的空间的必要性可打个比方来说。通常的认识是:人类生活在三维空间里。曾有传说,某人在某地方突然消失得踪迹全无,很快又在数千公里之外的某地出现。还有“特异功能者”能够从密封的药瓶里取出药片(而不必开启瓶盖或打破药瓶)。研究过多维空间的学者推理说:如果这是真的,那人和药片一定是利用了第四维空间的通道。由此看来,对于三维空间的凡人来说,视野及行动都能扩展至四维空间的,简直就是神仙了!同理,假设有一只生活在一维空间的小虫,它只能在一维的实数数轴上爬来爬去,它那可怜的视力根本看不到左右两侧;另一个二维空间的虫子却能在广阔的二维复平面上自由活动。 相对于一维的可怜虫,二维的小虫就是神仙。 许多经常做 FFT 分析计算的人往往忘记、或根本不知, FFT 程序的输入数组除了[ X ]之外,还有另一个[ Y ](尤其安装在一些大型工具软件系统里的 FFT ,当使用者把实信号数据输入[ X ]后,[ Y ]会自动设置为零)。国内有研究“故障诊断”的知名学者,在检测转轴 x 和 y 两垂直方向的径向误差后,将两组数据分别做 FFT (本来只需做一次 FFT ,却做了两次),再把得到的两个单边频谱合成为椭圆单边谱,并美其名曰“ 全息谱 ”或“ 全矢谱 ”。此法绕了弯路,数学推导并无错误,但是如此 舍近求远 的所谓创“新”却令人可叹又可笑。岂不知, Fourier 老先生早就为你的 x 、 y 信号预备好一个非常漂亮的双边频谱啦! 如果把 Fourier 变换比喻成一个双轮的手推车,人们历来竟一直把一侧轮胎悬空,使它歪着,当作独轮车来推。如今正在使用 FFT 的朋友,您还当它是个独轮车吗? Fourier变换复指数展开式的几何意义
个人分类: 未分类|18324 次阅读|6 个评论
FFT变换用EXCEL实现
FQOO 2009-9-14 21:57
EXCEL做FFT变换 用EXCEL函数做16点的FFT变换 试用一下博客 EXCEL密码:82607203
个人分类: 科研点滴|14761 次阅读|5 个评论

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-5-22 05:46

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部