CN107180640B - 一种相位相关的高密度叠窗频谱计算方法 - Google Patents

Abstract

本发明公开了一种相位相关的高密度叠窗频谱计算方法，包括：进行初始化，并输入原始信号；采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算，得到信号的频谱。本发明采用了基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法，通过在进行傅里叶变换时随源信号窗口的移动将原傅里叶变换矩阵的前若干行循环移动到最底端来使相邻窗重复的部分可以省略计算，节省了计算量，并采用了计算量小的IIR滤波器对窗函数进行拟合，计算量更低；频谱计算得到的信号的频谱包含有模值谱、相位谱和瞬时相位差谱，在保证模值分量正确检出的同时，能测得相位和相位差，更精确和直观。本发明可广泛应用于信号处理领域。

Description

一种相位相关的高密度叠窗频谱计算方法

技术领域

本发明涉及信号处理领域，尤其是一种相位相关的高密度叠窗频谱计算方法。

背景技术

频谱分析方式有着多样性的特征，就现阶段来看，频谱分析方式是多种多样的，有cross算法、DFT算法、prony算法、小波变换、卡尔曼滤波算法等，但是这些算法都存在各种各样的不足，其中，DFT算法的应用范围是最为广泛的，在高次谐波以及非整次谐波含量较少的情况下，该种算法的精度是十分理想的，该种算法应用了循环与递归算法，计算速度快，抗干扰性强，能够消除整次谐波对分析过程的不良影响。而使用加窗法与滤波法也能够避免DFT算法出现插值方向错误的问题。

对于一般的叠窗DFT算法来说，其窗口间隔会被设定为窗口长度的几分之一，而且当叠窗太疏时相位差不好用，从而导致一些基于叠窗FFT的应用会只保留模值，相位和相位差信息会被丢弃。但在实际应用中，目前的这种叠窗DFT算法会对音频处理产生以下不利的影响：由于受运算量所限，叠窗的间隔一般不能设置得太短，这样就只能从各个离散，稀疏的时间点观察数据。虽然对于模值而言，从各个离散和稀疏的时间点观察数据影响并不大，因为邻近窗口的模值是相似的，所以其测定值基本能反映实际状况；但其对于相位信息的影响就大为不同了，其可能会因设定了一个较长的窗口间隔而导致两个窗口之间的相位差超过π，这样就对相位信息变化速率的测定造成阻碍。

综上所述，业内亟需一种新的叠窗频谱计算方法，能降低运算量，并在保证模值分量正确检出的同时测得更精确和直观的相位变化信息。

发明内容

为解决上述技术问题，本发明的目的在于：提供一种计算量低、精确和直观的，相位相关的高密度叠窗频谱计算方法。

本发明所采取的技术方案是：

一种相位相关的高密度叠窗频谱计算方法，包括以下步骤：

进行初始化，并输入原始信号；

采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算，得到信号的频谱，所述基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法在进行傅里叶变换时随源信号窗口的移动将原傅里叶变换矩阵的前若干行循环移动到最底端来省略计算相邻窗的重复部分，并在加窗时采用了IIR滤波器对窗函数进行拟合，以节省计算量；其中，信号的频谱包含有模值谱、相位谱和瞬时相位差谱。

进一步，所述采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算，得到信号的频谱这一步骤，其包括：

Step1、使用原傅里叶变换矩阵M对原始信号S进行循环滚动卷积运算，得到中间计算矩阵并按照设定的存储格式存储中间计算矩阵；

Step2、根据窗函数w计算IIR参数k，并根据k对中间计算矩阵进行IIR滤波计算，得到近似拟合后的结果；

Step3、根据近似拟合后的结果得到相应信号的频谱。

进一步，所述步骤Step1包括：

Step11、定义原傅里叶变换矩阵M的前k行循环移动到最底端，得到移动后的傅里叶变换矩阵为M^(k)，其中，k为正整数；

Step12、定义S^(k)为抽取S的第k到k+n-1个元素所组成的一个行向量，记S^(k)＝(s_k,s_k+1,…,s_k+n-1)；

Step13、初始时，令k＝0；

Step14、对S^(k)与M^(k)进行卷积运算，并将卷积运算的中间结果按照设定的存储格式存储到中间计算矩阵中；

Step15、记窗口间隔为v，将M^(k)的前v行循环移动到M^(k)的最底端，得到M^(k+v)；

Step16、令k:＝k+v，使得步骤Step15中所述的M^(k+v)记为M^(k)，其中，“:＝”为赋值符号；

Step17、判断k+n-1是否已到达或超过S的长度，若是，则完成中间计算矩阵的计算过程；反之，则返回步骤Step14。

进一步，所述中间计算矩阵的存储规律为：中间计算矩阵中s同一列的下标相同，且s的下标从左到右各列递增1，中间计算矩阵中m同一行的上标相同，且m的上标从第一行到第n行逐行递增1，m同一列的下标相同，m的下标等于s下标模n，具体存储格式如下表所示：

上表中，S(s₀,s₁,…,s_i)为原始信号，s_i为S中第i个元素，i∈[0,t)，且i为自然数，记

其中，

表示M^(u)的第j列第i行元素，％为求余数运算，M^(u)为将原傅里叶变换矩阵M的前u行移动到最底端后得到的变换矩阵。

进一步，所述步骤Step2包括：

Step21、确定窗函数w；

Step22、根据窗函数w计算IIR参数k的值k₁,k₂,…,k_q；

Step23、对中间计算矩阵的每一行使用IIR参数k的值进行IIR滤波计算，得到结果矩阵A^*。

进一步，所述步骤Step23对中间计算矩阵的每一行能使用不同窗函数w计算出的不同IIR参数k来进行IIR滤波计算。

进一步，以中间计算矩阵的每一行作为处理对象，若将中间计算矩阵某行的第i个元素记为x_i，则步骤Step23包括：

Step231、根据IIR参数k的值对x_i按i＝0,i＝1,…,i＝t-1的正向顺序计算出x′_i，所述x′_i的计算公式为：

其中，q为IIR滤波器的阶数，且q远小于n，max为求最大值函数；

Step232、根据IIR参数k的值和x′_i按i＝t-1,i＝t-2,…,i＝0的反向顺序计算出x″_i，所述x″_i的计算公式为：

其中，x″_i为IIR滤波的近似拟合结果，min为求最小值函数；

Step233、将对应行的各x″_i的值写回中间计算矩阵x_i的相应位置并进行存储，最终得到结果矩阵A^*。

进一步，所述步骤Step3包括：

Step31、根据结果矩阵A^*计算相应信号的模值

以及相位

记

为A^*的第j行第i个元素，则模值

由

和

决定，

也由

和

决定，所述

和

的表达式分别为：

Step32、使用相位计算最终的相位差值D。

进一步，所述步骤Step32包括：

Step321、计算

中任意两个相邻元素的相位差值，所述

中任意两个相邻元素的相位差值

的计算公式为：

i∈[0,t-2]；

Step322、根据公式

进行逐项递推，计算出最终的相位差值D，其中，

为从下标i+1开始，横向相邻的r个相位差值的累加结果，i∈[0,t-1-r],r＜t，fixRot(x)为角度修正函数，fixRot(x)的表达式为：

进一步，所述采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算时，能只选择感兴趣频段进行频谱计算而省略不感兴趣频段的计算。

本发明的有益效果是：包括输入原始信号和采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算，得到信号的频谱的步骤，采用了基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法，通过在进行傅里叶变换时随源信号窗口的移动将原傅里叶变换矩阵的前若干行循环移动到最底端来使相邻窗重复的部分可以省略计算，节省了计算量，并采用了计算量小的IIR滤波器对窗函数进行拟合，计算量更低；频谱计算得到的信号的频谱包含有模值谱、相位谱和瞬时相位差谱，在保证模值分量正确检出的同时，能测得相位和相位差，更精确和直观。进一步，对原始信号进行频谱计算时，能只选择感兴趣频段进行频谱计算而省略不感兴趣频段的计算，提高了计算效率。

附图说明

图1为本发明一种相位相关的高密度叠窗频谱计算方法的具体流程图；

图2为传统傅里叶变换方法的运算形式示意图；

图3为实施例一高密度叠窗算法矩阵的移动过程示意图；

图4为实施例一高密度叠窗算法矩阵的预乘优化过程示意图；

图5为实施例一采用本发明的高密度叠窗算法计算低频部分时的运算过程示意图。

具体实施方式

参照图1，一种相位相关的高密度叠窗频谱计算方法，包括以下步骤：

进行初始化，并输入原始信号；

采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算，得到信号的频谱，所述基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法在进行傅里叶变换时随源信号窗口的移动将原傅里叶变换矩阵的前若干行循环移动到最底端来省略计算相邻窗的重复部分，并在加窗时采用了IIR滤波器对窗函数进行拟合，以节省计算量；其中，信号的频谱包含有模值谱、相位谱和瞬时相位差谱。

本发明在进行变换矩阵循环滚动和IIR滤波器拟合后，还需要进行一些额外的运算(如图1所示的相位差计算等)才能得到最终的模值、相位和瞬时相位差。

参照图1、3和4，进一步作为优选的实施方式，所述采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算，得到信号的频谱这一步骤，其包括：

Step1、使用原傅里叶变换矩阵M对原始信号S进行循环滚动卷积运算，得到中间计算矩阵并按照设定的存储格式存储中间计算矩阵；

Step2、根据窗函数w计算IIR参数k，并根据k对中间计算矩阵进行IIR滤波计算，得到近似拟合后的结果；

Step3、根据近似拟合后的结果得到相应信号的频谱。

参照图3和4，进一步作为优选的实施方式，所述步骤Step1包括：

Step11、定义原傅里叶变换矩阵M的前k行循环移动到最底端，得到移动后的傅里叶变换矩阵为M^(k)，其中，k为正整数；

Step12、定义S^(k)为抽取S的第k到k+n-1个元素所组成的一个行向量，记S^(k)＝(s_k,s_k+1,…,s_k+n-1)；

Step13、初始时，令k＝0；

Step14、对S^(k)与M^(k)进行卷积运算，并将卷积运算的中间结果按照设定的存储格式存储到中间计算矩阵中；

Step15、记窗口间隔为v，将M^(k)的前v行循环移动到M^(k)的最底端，得到M^(k+v)；

Step16、令k:＝k+v，使得步骤Step15中所述的M^(k+v)记为M^(k)，其中，“:＝”为赋值符号；

Step17、判断k+n-1是否已到达或超过S的长度，若是，则完成中间计算矩阵的计算过程；反之，则返回步骤Step14。

其中，各种循环滚动操作是逻辑上的操作，在实际计算中，M到M^(k)的变换仅仅通过改变对M的读取顺序来实现，并不实际改变M的存储值。

进一步作为优选的实施方式，所述中间计算矩阵的存储规律为：中间计算矩阵中s同一列的下标相同，且s的下标从左到右各列递增1，中间计算矩阵中m同一行的上标相同，且m的上标从第一行到第n行逐行递增1，m同一列的下标相同，m的下标等于s下标模n，具体存储格式如下表所示：

上表中，S(s₀,s₁,…,s_i)为原始信号，s_i为S中第i个元素，i∈[0,t)，且i为自然数，记

其中，

表示M^(u)的第j列第i行元素，％为求余数运算，M^(u)为将原傅里叶变换矩阵M的前u行移动到最底端后得到的变换矩阵。

参照图1，进一步作为优选的实施方式，所述步骤Step2包括：

Step21、确定窗函数w；

Step22、根据窗函数w计算IIR参数k的值k₁,k₂,…,k_q；

Step23、对中间计算矩阵的每一行使用IIR参数k的值进行IIR滤波计算，得到结果矩阵A^*。

本发明可采用Butterworth，或Chebyshev等现有的低通滤波器设计方法来根据窗函数w确定IIR参数k的值，也就是说IIR参数k的值和w的转换关系与IIR滤波器的具体设计有关。

进一步作为优选的实施方式，所述步骤Step23对中间计算矩阵的每一行能使用不同窗函数w计算出的不同IIR参数k来进行IIR滤波计算。

其中，步骤Step23中间计算矩阵的每一行，可使用不同的w，以确定不同的IIR参数k，以模拟使用不同σ参数的高斯窗或其它类型宽窄不同的窗，从而让不同频段拥有不同的时间响应特性。

进一步作为优选的实施方式，以中间计算矩阵的每一行作为处理对象，若将中间计算矩阵某行的第i个元素记为x_i，则步骤Step23包括：

Step231、根据IIR参数k的值对x_i按i＝0,i＝1,…,i＝t-1的正向顺序计算出x′_i，所述x′_i的计算公式为：

其中，q为IIR滤波器的阶数，且q远小于n，max为求最大值函数；

Step232、根据IIR参数k的值和x′_i按i＝t-1,i＝t-2,…,i＝0的反向顺序计算出x″_i，所述x″_i的计算公式为：

其中，x″_i为IIR滤波的近似拟合结果，min为求最小值函数；

Step233、将对应行的各x″_i的值写回中间计算矩阵x_i的相应位置并进行存储，最终得到结果矩阵A^*。

参照图1，进一步作为优选的实施方式，所述步骤Step3包括：

Step31、根据结果矩阵A^*计算相应信号的模值

以及相位

记

为A^*的第j行第i个元素，则模值

由

和

决定，

也由

和

决定，所述

和

的表达式分别为：

Step32、使用相位计算最终的相位差值D。

进一步作为优选的实施方式，所述步骤Step32包括：

Step321、计算

中任意两个相邻元素的相位差值，所述

中任意两个相邻元素的相位差值

的计算公式为：

i∈[0,t-2]；

Step322、根据公式

进行逐项递推，计算出最终的相位差值D，其中，

为从下标i+1开始，横向相邻的r个相位差值的累加结果，i∈[0,t-1-r],r＜t，fixRot(x)为角度修正函数，fixRot(x)的表达式为：

进一步作为优选的实施方式，所述采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算时，能只选择感兴趣频段进行频谱计算而省略不感兴趣频段的计算。

下面结合说明书附图和具体实施例对本发明作进一步解释和说明。

实施例一

现有的叠窗DFT算法当叠窗太疏时相位差往往会因不好用而被丢弃，会对相位信息变化速率的测定造成阻碍；而当叠窗密度高时，其运算量又会很大。为了解决上述问题，本发明提出了一种相位相关的高密度叠窗频谱计算方法，通过改变傅里叶变换的传统计算方法，使其能更好地对数据进行高密度傅里叶叠窗计算，在保证模值分量准确算出的同时，测得更精确和直观的瞬时相位变化信息。

本发明相位相关的高密度叠窗频谱计算方法的具体实现过程为：

(一)进行初始化，并输入原始信号。

其中，原始信号根据实际需要进行调整，其主要包括PCM信号和音频处理信号等信号。

如图1所示，此过程可进一步细分为如下步骤：

(1)确定变换矩阵的维度n*n；

(2)构建原傅里叶变换矩阵M；

(3)输入原始PCM信号S。

输入原始信号后，接着可采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算。

(二)PFT算法原理

为了克服传统叠窗DFT算法的缺陷，本实施例提出了基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法(简称PFT算法)。

PFT算法的核心计算方式与原始的傅里叶变换(即FT)近似，通过对原始信号乘以不同频率的sin和cos函数从而得到信号不同频率的分量。但PFT算法与FT的有以下几个不同点：

1)sin和cos函数的初始相位随着窗口的移动而改变，体现为矩阵M的循环移动；

2)为抑制叠窗FT的栅栏效应，FT的一般做法是对信号先进行乘窗函数处理后再进行傅里叶变换，而PFT算法的加窗过程与傅里叶变换是同时进行的；

3)PFT算法在运算过程中，相近的窗之间大部分的运算过程和运算结果是共享的，每一个窗的平均运算量很小。

(三)快速叠窗过程

(1)快速叠窗过程概述

传统的FT运算形式如图2所示，所以实际上FT就是将原始数据S分窗口乘以变换矩阵M得到中间结果(a₀,b₀,a₁,b₁,a₂,b₂…a_n/2-1,b_n/2-1)^T，n为变换矩阵M的阶数。为更清楚地说明S分窗口的细节，记S^(k)为抽取S的第k到k+n-1个元素所组成的一个行向量，即有S^(k)＝(s_k,s_k+1,…,s_k+n-1)；

定义计数变量k，用于记录当前计算的窗口的最左端元素对应S的下标。初始时，k＝0；记PFT的窗口间隔为v，n为变换矩阵M的阶数。

对S^(k)与M^(k)进行点积运算，得到一个窗口的中间结果，记为列向量A^(k/v)＝(a₀,b₀,a₁,b₁,a₂,b₂…a_n/2-1,b_n/2-1)^T。

不同的k得到的各个中间结果依次横向并行排列，组合成中间计算矩阵A。

如果直接用FT来做叠窗变换，由于不同的窗的S不同，所以每个窗都要做图2的计算，这样运算量会很大；而且当叠窗间隔很短时，即使采用FFT来做运算，其运算量也会相当大，因此需要对此做优化。对图2的M进行改造，每计算完一次S^(k)与M^(k)的点积后，下一个窗口将为S^(k+v)，若在计算窗口S^(k+v)时，将M^(k)的前v行循环移动到最底端，将移动后的结果记为M^(k+v)，如图3所示，则A^(k)＝S^(k)×M^(k)与A^(k+v)＝S^(k+v)×M^(k+v)中的A^(k)和A^(k+v)存在较大比例的完全一致的部分，如图4的斜线阴影部分所示；而不同的需要计算的部分，如图4的竖线阴影部分所示，则只占少数。具体而言，有

i∈(0,n-u)，因此，对于

而言，下标在i∈(0,n-u)内的部分直接使用

的结果即可，不需要进行计算。所以，窗口间隔v越短，A^(k)和A^(k+v)中相同的区域就越大，相对而言，每个窗口所需要的计算量就越小。

从使用窗函数的角度出发，应对中间计算矩阵A的每一行进行高斯模糊，并拟合高斯模糊的近似结果覆盖原来的矩阵A。然而，高斯模糊是卷积操作，计算量很大，所以改为使用IIR滤波器来近似拟合高斯模糊的结果，使用IIR滤波后的结果覆盖A原来的数据，并将IIR滤波后的结果记为A^*。实际上这两种计算方法结果很相似，最大的不同在于计算量。所以后文中在论及计算结果的区别时并不严格区分此步骤具体使用了这两种方法中的何种方法。

若前面按传统的FT计算方法，不对M进行滚动操作而直接做卷积，然后得到高斯模糊的结果，并用得到高斯模糊的结果覆盖原来的矩阵A来进行近似拟合，将高斯模糊后的A记为A^#，将A^#转换为相应的模值谱和相位谱，则有：模值谱

相位谱

若前面按A^(k)＝S^(k)×M^(k)进行计算，将A做IIR滤波得到A^*，使用A^*求模值谱和相位谱，则模值谱

相位谱

上面两种计算方式所得到的相位谱Λ，模值谱Γ和相位谱Λ^*，模值谱Γ^*，有相同的物理意义。故相位谱Λ，模值谱Γ和相位谱Λ^*，模值谱Γ^*物理意义相同，其中，Γ≡Γ^*，而Λ^*则与Λ中的每一个对应的分量都相差一个确定的角度，这个角度差仅与下标i相关。虽然A^*与A不相等，但是两个结果可以互相转换；且实际的应用测试表明，往往A^*，Γ^*，Λ^*的结果直接用起来更方便。

综上，使用高斯模糊来替换加窗操作，是一种等效的变换方法，本身并不改变运算量，但有利于后面使用简化算法。简化算法分别使用了A^(k)＝S^(k)×M^(k)代替传统FT和使用了IIR滤波器替代高斯模糊，构成了快速计算高密度叠窗FT方法的两个方面，这两个方法各自都能提高计算效率，且可以组合使用。值得一提的是，在上述的运算过程中，不需要每一步去改写M，只需要在每次做乘法的时候改变对M的读取顺序即可。

(2)快速叠窗过程扩展

如图4所示，斜纹阴影部分用于计算A^(u)和A^(u+v)中乘法计算重复的部分。可见当v较短时，每个窗口都可以从上一个窗口中得到绝大部分的中间结果，从而在计算A^(u+v)时只需要更新图4竖纹阴影部分所示的少量乘积，就可以得到本窗口的变换结果。

此外，若只对部分频段感兴趣，例如只对低频部分感兴趣，则可以只计算低频部分的数据。如图5所示，图5只要计算斜纹阴影部分的乘积，就可以得到竖纹部分的结果，这在一些应用场合，例如语音信号处理中可以使用。

(3)快速叠窗过程的计算推理细节

在实际的时域-频域变换处理过程中，除了上面提到的做乘法部分运算外，还要做乘以窗函数和将各个乘积相加的工作。为了使这部分的运算更简便，本实施例对上述步骤的中间乘法结果的存储方式进行了如下表1的约定：

表1中间计算矩阵

而当v为1时，无窗函数的变换处理结果矩阵如表2所示：

表2无窗函数的变换处理结果矩阵

由表2可知，若要直接计算所有的求和操作，其运算量将会很大，但是，由于

即式(1)左边的和式可由右边的公式推出，所以实际的计算量将会降低很多。

当窗口间隔v为p时，递推式与上式(1)近似，同样可以有效减少运算量：

为了得到更清晰的结果，在进行叠窗傅里叶变换时，常会通过预乘窗函数来抑制“栅栏效应”。而在v为1时，加窗函数的变换处理结果矩阵如下表3所示。

表3加窗函数的变换处理结果矩阵

表3中，w表示窗函数，但由于有窗函数w的存在，求和公式不能再像公式(1)和公式(2)那样进行递推了。虽然w具体是什么形状会对旁瓣大小和主峰精确度等造成一定影响，但只要w是一个中间高两边低的函数，就可达到减少栅栏效应的目的。所以w有很大的选择空间。利用这一特性，可以选取合适的w，使表3可以进行快速近似计算。

不失一般性，此处w选取正态分布函数

来作为窗函数，其中，均值

标准差σ≤n/6，

为向下取整符号。

从高斯模糊的概念而言，对信号做高斯模糊就是使用式(3)的正态分布函数对信号进行卷积，而正态分布有3σ特性，即

所以，虽然w(x)的定义域是无穷的，但是在实际做卷积时，只需要使用中心的±3σ的长度进行卷积就可达到非常接近的效果。而由于前面选择的σ≤n/6，所以，对于表3而言，如果将某一行抽出来看，可以发现其结果近似相当于对表2对应的行进行了高斯模糊。高斯模糊本质上是一个卷积过程，如果不对正态分布函数进行合理截断，其运算量甚至将远远超过表3中的运算量；即使进行了截断，其运算量也不会比表3的运算量有明显改善，所以直接使用高斯模糊并不可行。

但是，可以参照这个思路，使用计算量很小的IIR滤波器来进行拟合，这是最关键的一步。具体而言，将中间矩阵某行的第i个元素抽出，记为x_i，i∈[0,t)，如果做传统的高斯卷积就是

然而这一过程的计算量比较大，为减少运算量，可以构造合理的k_j取值，按i＝0,i＝1,…,i＝t-1的正向顺序先执行如下的计算过程：

再按i＝t-1,i＝t-2…,i＝0的反向顺序执行如下的计算过程：

从而使算得的近似结果

以上就是IIR滤波器的思路。在k_j设置合理的前提下，q越大，

和x″_i的拟合度就越高。

对于IIR参数k，可采用但不限于Butterworth，或Chebyshev等现有的低通滤波器设计方法来确定，也就是说IIR参数k的值和w的转换关系与IIR滤波器的具体设计有关。

值得注意的是，使用IIR滤波器时，q一般取很小的值，一般有q远小于n。在处理长窗口时，例如窗口超过512时，低阶的算法会产生抖动问题，所以在实际使用长窗口时，需要对其进行改进，适当使用更高阶数的IIR滤波器。但是即使窗口再长，在处理音频信号时，其阶数也一般不会超过10阶。所以，在使用IIR滤波器的近似算法对卷积过程进行拟合后，其运算量将明显小于直接进行高斯模糊的运算量。

在实际使用时，用于计算x″_i的x_i其实就是表2中的

为了后面说明的方便，在此增加行标j，将第j行的x″_i记作

本发明通过两个元素通过求平方和再开方的方式得到模值，通过这两个元素的反正切函数来得到相位。从第0行开始，每两行中的一对

值确定一对模值和相位，则有：

实际上，前面高斯模糊所使用的正态分布窗只是一个特例，若要使用不同形状的窗口来获得不同的旁瓣主峰特性，则可以通过改变窗函数w来调整IIR滤波器的参数k的方式来实现。

IIR滤波器计算x’和x”时是可以在原位进行结果替换的，所以IIR滤波器相比起传统高斯模糊而言，空间的要求也更小。

另外，由于IIR滤波器是无限响应的，所以实际上窗口是无限长的，但是可以通过调节窗函数主峰的宽度来获得等效窗口长度。例如，在使用IIR滤波器模拟正态窗函数时，由于±3σ中已包含了99％以上的面积，所以可以把窗口长度看作是±3σ，在这个情况下，可以把±3σ称作等效窗口长度。所以实际上等效窗口长度是由窗函数来决定的。

由于是做横向模糊，实际上对于不同的频率(对应表1中的不同行)，可以采用不同的IIR参数，以模拟使用不同σ参数的高斯窗或其它类型宽窄不同的窗，从而让不同频段拥有不同的时间响应特性。例如，对信号的低频分量使用较长等效窗口长度(对应高斯窗的σ较大)的窗函数，能保证低频分量能被很好地检出；而对高频分量施以较短等效窗口长度(对应高斯窗的σ较小)的窗函数，在保证高频分量检出的同时又能更清晰地检出高频分量的瞬时变化。从这一点看，PFT算法比普通叠窗傅里叶变换要更加灵活。

(四)计算信号的相位差。

确定A近似拟合后的结果矩阵A^*后，即可得到A^*的模值

以及相位

接着，可根据Λ^*进行相位差计算，具体计算过程为：

记角度修正函数

定义相位差为

和

所在的时间点相隔很近，且自然声音，尤其是语音在瞬时内相位差不会有明显变化，所以，为了得到更加精确的相位差，可以计算短时间段内的相位差和值作为计算结果(即为取得更平滑结果，可以取横向相邻的r个相位差结果做累加)，则可根据下式计算最终的相位差

其中，i∈[0,t-1-r],r＜t

从上式可见：

所以在实际计算时，可以预先计算所有

i∈[0,t-2]，然后通过

快速计算最终的相位差D。

通过本发明的方法计算出信号的模值、相位和相位差，即可根据信号的模值、相位和相位差进行后续的频谱分析。

本发明提出了一种相位相关的高密度叠窗频谱计算方法，通过改变傅里叶变换的传统计算方法，使其能更好地对数据进行高密度傅里叶叠窗计算，在保证模值分量准确算出的同时，测得更精确和直观的瞬时相位变化信息。与现有技术相比，本发明的方法具有以下优点：

(1)能测得相位和幅值的精确连续变化情况；

(2)能有目的地选择频段，对不关心的频段可以省略计算；

(3)对不同频段可以用不同的窗口间隔来分析，典型情况是：对低频信号使用较长的窗口，而对高频信号则使用较短的窗口，这样既可以准确识别低频分量，又能对高频信号获得较好的时间响应特性；

(4)对不同频段可以使用不同的频率间隔，例如在处理语音识别信号时，可以对高频部分施以更宽的频率间隔，这既符合人类的听觉特性，又能节省计算量。

以上是对本发明的较佳实施进行了具体说明，但本发明并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可做作出种种的等同变形或替换，这些等同的变形或替换均包含在本申请权利要求所限定的范围内。

Claims (10)

1.一种相位相关的高密度叠窗频谱计算方法，其特征在于：包括以下步骤：

进行初始化，并输入原始信号；

采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算，得到信号的频谱，所述基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法在进行傅里叶变换时随源信号窗口的移动将原傅里叶变换矩阵的前若干行循环移动到最底端来省略计算相邻窗的重复部分，并在加窗时采用了IIR滤波器对窗函数进行拟合，以节省计算量；其中，信号的频谱包含有模值谱、相位谱和瞬时相位差谱。

2.根据权利要求1所述的一种相位相关的高密度叠窗频谱计算方法，其特征在于：所述采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算，得到信号的频谱这一步骤，其包括：

Step1、使用原傅里叶变换矩阵M对原始信号S进行循环滚动卷积运算，得到中间计算矩阵并按照设定的存储格式存储中间计算矩阵；

Step2、根据窗函数w计算IIR参数k，并根据k对中间计算矩阵进行IIR滤波计算，得到近似拟合后的结果；

Step3、根据近似拟合后的结果得到相应信号的频谱。

3.根据权利要求2所述的一种相位相关的高密度叠窗频谱计算方法，其特征在于：所述步骤Step1包括：

Step11、定义原傅里叶变换矩阵M的前k行循环移动到最底端，得到移动后的傅里叶变换矩阵为M^(k)，其中，k为正整数；

Step12、定义S^(k)为抽取S的第k到k+n-1个元素所组成的一个行向量，记S^(k)＝(s_k,s_k+1,…,s_k+n-1)；

Step13、初始时，令k＝0；

Step14、对S^(k)与M^(k)进行卷积运算，并将卷积运算的中间结果按照设定的存储格式存储到中间计算矩阵中；

Step15、记窗口间隔为v，将M^(k)的前v行循环移动到M^(k)的最底端，得到M^(k+v)；

Step16、令k:＝k+v，使得步骤Step15中所述的M^(k+v)记为M^(k)，其中，“:＝”为赋值符号；

Step17、判断k+n-1是否已到达或超过S的长度，若是，则完成中间计算矩阵的计算过程；反之，则返回步骤Step14。

4.根据权利要求3所述的一种相位相关的高密度叠窗频谱计算方法，其特征在于：所述中间计算矩阵的存储规律为：中间计算矩阵中s同一列的下标相同，且s的下标从左到右各列递增1，中间计算矩阵中m同一行的上标相同，且m的上标从第一行到第n行逐行递增1，m同一列的下标相同，m的下标等于s下标模n，具体存储格式如下表所示：

上表中，S(s₀,s₁,…,s_i)为原始信号，s_i为S中第i个元素，i∈[0,t)，且i为自然数，记

其中，

表示M^(u)的第j列第i行元素，％为求余数运算，M^(u)为将原傅里叶变换矩阵M的前u行移动到最底端后得到的变换矩阵。

5.根据权利要求4所述的一种相位相关的高密度叠窗频谱计算方法，其特征在于：所述步骤Step2包括：

Step21、确定窗函数w；

Step22、根据窗函数w计算IIR参数k的值k₁,k₂,…,k_q；

Step23、对中间计算矩阵的每一行使用IIR参数k的值进行IIR滤波计算，得到结果矩阵A^*。

6.根据权利要求5所述的一种相位相关的高密度叠窗频谱计算方法，其特征在于：所述步骤Step23对中间计算矩阵的每一行能使用不同窗函数w计算出的不同IIR参数k来进行IIR滤波计算。

7.根据权利要求6所述的一种相位相关的高密度叠窗频谱计算方法，其特征在于：以中间计算矩阵的每一行作为处理对象，若将中间计算矩阵某行的第i个元素记为x_i，则步骤Step23包括：

Step231、根据IIR参数k的值对x_i按i＝0,i＝1,…,i＝t-1的正向顺序计算出x′_i，所述x′_i的计算公式为：

其中，q为IIR滤波器的阶数，且q远小于n，max为求最大值函数；

Step232、根据IIR参数k的值和x′_i按i＝t-1,i＝t-2,…,i＝0的反向顺序计算出x″_i，所述x″_i的计算公式为：

其中，x″_i为IIR滤波的近似拟合结果，min为求最小值函数；

Step233、将对应行的各x″_i的值写回中间计算矩阵x_i的相应位置并进行存储，最终得到结果矩阵A^*。

8.根据权利要求6所述的一种相位相关的高密度叠窗频谱计算方法，其特征在于：所述步骤Step3包括：

Step31、根据结果矩阵A^*计算相应信号的模值

以及相位

记

为A^*的第j行第i个元素，则模值

由

和

决定，

也由

和

决定，所述

和

的表达式分别为：

Step32、使用相位计算最终的相位差值D。

9.根据权利要求8所述的一种相位相关的高密度叠窗频谱计算方法，其特征在于：所述步骤Step32包括：

Step321、计算

中任意两个相邻元素的相位差值，所述

中任意两个相邻元素的相位差值

的计算公式为：

Step322、根据公式

进行逐项递推，计算出最终的相位差值D，其中，

为从下标i+1开始，横向相邻的r个相位差值的累加结果，i∈[0,t-1-r],r＜t，fixRot(x)为角度修正函数，fixRot(x)的表达式为：

10.根据权利要求1-9任一项所述的一种相位相关的高密度叠窗频谱计算方法，其特征在于：所述采用基于变换矩阵循环滚动和IIR滤波器拟合的高密度叠窗算法对原始信号进行频谱计算时，能只选择感兴趣频段进行频谱计算而省略不感兴趣频段的计算。

Images