CN104050147A - 将时域信号转换成频域信号的方法与系统 - Google Patents

Abstract

本发明涉及数字信号处理领域，公开了一种将时域信号转换成频域信号的方法与系统。该方法将时域信号用曲线平滑后再生为近似模拟信号，传递给进行有阻尼的受迫简谐振动的谐振系统，并计算该谐振系统的谐振强度，以作为频域信号强度。因为可以任选频点计算，所以在只要计算特定的一个或几个频点的信号时，计算量远低于传统计算方法。进一步通过计算频域信号强度变化过程，得到时间-频率分布图，若将相应频点的所有采样点对应时刻的频域信号强度叠加即可得到该频点的频域信号强度，从而成功将时域信号转换成频域信号。由于是逐个采样点计算，时间响应最高可达到采样周期并且响应幅度可调。

Description

将时域信号转换成频域信号的方法与系统

技术领域

本发明涉及数字信号处理领域，特别涉及将时域信号转换成频域信号的方法与系统。

背景技术

快速傅里叶变换是目前对数字信号数据最常用的处理方法之一，它可以将信号从时域变换到频域，产生统计频谱，进行进一步研究。而为了对诸如语音信息等频谱不断变化的信号进行研究，发展了时频变换方法。时频变换方法目前一些主要的方法是短时傅里叶变换和小波变换，以及魏格纳-维尔变换等，但最主要的方法是短时傅里叶变换。

短时傅里叶变换(Short Time Fourier Transform)具体做法是将数字信号划分成固定长度的帧，相邻帧之间有50%～90%的重叠，对每一帧数据先加窗处理，然后做快速傅里叶变换(Fast Fourier Transform)得到该帧数据的一个频率域分布，将所有帧的频率域分布按时间顺序排列，就得到了三维的时间-频率分布图。

短时傅里叶变换具有速度快，计算量小的特点，但是其产生的时间-频率分布图，频率的分割点数永远等于数据帧长度的一半，是均匀分布的，并且由于功率值反映的是每帧数据对应时间片内的平均值，对瞬时变化的时间响应精度只能达到每帧时长的1/2-1/5。因此，给定一个采样频率后，变换所得到时频分布图的时间精度和频率精度的乘积基本为定值，我们只能在时间精度和频率精度之间被迫作折衷选择，难以得到同时具有高频率精度和高时间精度的时-频分布图。

对人耳蜗的研究表明，耳蜗有大约24000个听觉传感器，每个传感器都对特定频率的信号产生强烈反应并将信号幅度转换成电信号，众多听觉传感器通过神经束以并行的方式将这些电信号传递到大脑，大脑处理这些信号产生听觉。所以，从理论上我们也应该可以有一种类似听觉传感器的方法对信号进行并行的时频变换，能够在时间响应和频率精度上均高于短时傅里叶变换的转换方法。

发明内容

本发明的目的在于提供一种将时域信号转换成频域信号的方法与系统，时间响应最高可达到采样周期并且响应幅度可调。

为解决上述技术问题，本发明的实施方式公开了一种将时域信号转换成频域信号的方法，包括以下步骤：

根据时域信号中前一采样点和当前采样点确定平滑曲线f(t)的函数形式和相应参数；

将f(t)的参数和前一采样点对应时刻的状态参数输入到一个偏微分方程的解析解中，求得当前采样点对应时刻的状态参数，其中该偏微分方程代表一个有阻尼的受迫谐振系统，形式为k是该谐振系统无阻尼时的谐振角频率，n为该谐振系统的阻尼因子，k和n对应该谐振系统的一个共振角频率ω_r，状态参数中包含该谐振系统的速度v和位移x；

根据当前采样点对应时刻的状态参数中的速度v和位移x计算谐振强度，作为频域中与共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度。

本发明的实施方式还公开了一种声纹识别方法，在声纹识别的过程中使用上文的方法根据声音信号产生时间-频率分布图。

本发明的实施方式还公开了一种地震波分析方法，使用上文的方法根据地震波产生时间-频率分布图，以对地震波进行时频分析。

本发明的实施方式还公开了一种电力线路谐波分析方法，使用上文的方法根据电力线路谐波产生时间-频率分布图，以对电力线路谐波进行分析。

本发明的实施方式还公开了一种雷达的回波分析方法，使用上文的方法根据雷达的回波信号产生时间-频率分布图，以对回波信号进行时频分析。

本发明的实施方式还公开了一种将时域信号转换成频域信号的系统，包括：

平滑模块，用于根据时域信号中前一采样点和当前采样点确定平滑曲线f(t)的函数形式和相应参数；

谐振模块，用于将f(t)的参数和前一采样点对应时刻的状态参数输入到一个偏微分方程的解析解中，求得当前采样点对应时刻的状态参数，其中该偏微分方程代表一个有阻尼的受迫谐振系统，形式为 k是该谐振系统无阻尼时的谐振角频率，n为该谐振系统的阻尼因子，k和n对应该谐振系统的一个共振角频率ω_r，状态参数中包含该谐振系统的速度v和位移x；

计算模块，用于根据当前采样点对应时刻的状态参数中的速度v和位移x计算谐振强度，作为频域中与共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度。

本发明实施方式与现有技术相比，主要区别及其效果在于：

将时域信号用曲线平滑后再生为近似模拟信号，传递给进行有阻尼的受迫简谐振动的谐振系统，并计算该谐振系统的谐振强度，以作为频域信号强度，因为是采用特定形式曲线来作为外强迫力，从而可以直接采用解析解进行精确计算，避免了用诸如欧拉法、威尔逊-θ法或线性加速法等积分算法来求近似解，在提高计算准确性的同时极大地降低了计算量，并且由于可以任选频点计算，所以在只要计算特定的一个或几个频点的信号时，计算量远低于像快速傅里叶变换这类传统计算方法。

该系统包括平滑模块、谐振模块和计算模块，其中该平滑模块将时域信号用曲线平滑后再生为近似模拟信号，传递给该谐振模块进行有阻尼的受迫简谐振动，计算模块计算该受迫简谐振动的谐振强度，以作为频域信号强度，由于可以任选频点计算，所以在只要计算特定的一个或几个频点的信号时，计算量远低于传统计算方法，并且由于是可以直接采用解析解进行计算，从而不需要用诸如欧拉法、威尔逊-θ法或线性加速法等积分算法来求近似解，在保持计算准确性的同时极大地降低了计算复杂度。

进一步地，通过计算频域信号强度变化过程，进而得到时间-频率分布图，若将相应频点的所有采样点对应时刻的频域信号强度叠加即可得到该频点的频域信号强度，从而成功将时域信号转换成频域信号，因为可以任选频点计算，频率精度与采样频率无关，所以即使时域信号采样频率低，频率精度也可以远高于快速傅里叶变换，而且还可以任意选择感兴趣的频率区间来设置任意多个频率点进行转换，由于是逐个采样点计算，时间响应最高可接近采样周期，并且响应幅度可调，因此可以得到同时具有高频率精度和高时间精度的时间-频率分布图。

进一步地，对时域信号中各采样点或具有不同共振角频率的谐振系统采用相同或不同的曲线函数形式进行平滑连接，可适应一些特殊需求的场合。

进一步地，通过对共振角频率ω_r对应的受迫谐振系统的阻尼因子n进行调整，可以提高受迫谐振系统的频域分辨精度或时域分辨精度。

进一步地，通过不同方法得到与各共振角频率ω_r对应的谐振系统的阻尼因子n和谐振系统无阻尼时的谐振角频率k，可根据实际需求选择相应的谐振系统。

附图说明

图1是本发明第一实施方式中一种将时域信号转换成频域信号的方法的流程示意图；

图2是本发明第一实施方式中有阻尼的谐振系统在受到外强迫力持续作用后系统的能量变化的时间特性曲线；

图3是本发明第一实施方式中有阻尼的谐振系统在受到外强迫力持续作用时系统的谐振强度的频率特性曲线；

图4是本发明第一实施方式中采样值点由平滑曲线分段连接的示意图；

图5是本发明第二实施方式中一种将时域信号转换成频域信号的方法的流程示意图；

图6是本发明第三实施方式中一种将时域信号转换成频域信号的方法的流程示意图；

图7a至7g是本发明第四实施方式中对语音片段、粉红噪声、调频调幅混合声，分别由短时傅里叶变换和本发明转换方法进行转换得到的频谱图结果；

图8是本发明第八实施方式中一种将时域信号转换成频域信号的系统的结构示意图；

图9是本发明第九实施方式中一种将时域信号转换成频域信号的系统的结构示意图。

具体实施方式

在以下的叙述中，为了使读者更好地理解本申请而提出了许多技术细节。但是，本领域的普通技术人员可以理解，即使没有这些技术细节和基于以下各实施方式的种种变化和修改，也可以实现本申请各权利要求所要求保护的技术方案。

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明的实施方式作进一步地详细描述。

本发明第一实施方式涉及一种将时域信号转换成频域信号的方法。图1是该将时域信号转换成频域信号的方法的流程示意图。

具体地说，如图1所示，该将时域信号转换成频域信号的方法包括以下步骤：

在步骤101中，根据时域信号中前一采样点和当前采样点确定平滑曲线f(t)的函数形式和相应参数。

可以理解，f(t)可以对于时域信号中各采样点和具有不同共振角频率ω_r的谐振系统采用相同平滑曲线函数形式，也可以对于时域信号中各采样点或具有不同共振角频率ω_r的谐振系统采用不同平滑曲线函数形式。

例如，f(t)可以采用a+b cos(pt)的形式，其中p为时域信号的采样角频率的一半，a为前一采样值与当前采样值的算数平均值，b为前一采样值与当前采样值的差的一半，可以采用a′+b′p′t的形式，其中p′为时域信号的采样频率，a′为前一采样值，b′为前一采样值与当前采样值的差，可以采用a₀+a₁p′t+a₂p′²t²+a₃p′³t³的形式，其中p′为时域信号的采样频率，也可以采用a+bcos(pt)、a’+b’p’t和a₀+a₁p′t+a₂p′²t²+a₃p′³t³的形式组合。

此外，可以理解，在本发明的其他实施方式中，f(t)还可以采用其他的平滑曲线函数形式。

对时域信号中各采样点或具有不同共振角频率的谐振系统采用相同或不同的曲线函数形式进行平滑连接，可适应一些特殊需求的场合。

需注意的是，这些平滑曲线的参数可以通过直接输入前后两个或多个采样点的值来确定，也可以只将相关采样点输入到一个算法，由该算法进行预定曲线形式的曲线拟合或平滑，最终得到曲线参数，例如三次样条曲线、三次拟合曲线等。只要这些平滑曲线可以通过或接近波形图中前一采样点和当前采样点就可以了。

此后进入步骤102，将f(t)的参数和前一采样点对应时刻的状态参数输入到一个偏微分方程的解析解中，求得当前采样点对应时刻的状态参数。其中该偏微分方程代表一个有阻尼的受迫谐振系统，形式为 k是该谐振系统无阻尼时的谐振角频率，n为该谐振系统的阻尼因子，k和n对应该谐振系统的一个共振角频率ω_r，状态参数中包含该谐振系统的速度v和位移x。

可以理解，求一个偏微分方程的解析解是一种成熟的数字方法，在教科书中可以找到，本申请中不进行详细说明了。

此外，在本发明的各个实施方式中，可以通过以下几个方式来得到谐振系统的k和n：

1）谐振系统的阻尼因子n不变且为大于零的实数，由得到与各共振角频率ω_r对应的谐振系统无阻尼时的谐振角频率k；

2）谐振系统的振幅衰减倍数C不变且为大于零的实数，由得到与各共振角频率ω_r对应的谐振系统无阻尼时的谐振角频率k，再由得到与各共振角频率ω_r对应的阻尼因子n，其中振幅衰减倍数C为f(t)的振幅与谐振系统的谐振振幅A的比例，即定义稳态时表示外强迫力的f(t)的振幅与谐振系统模型因为共振而产生的谐振振幅A的比例定义为振幅衰减倍数C；

3）谐振系统的能量衰减倍数Q不变且为大于零的实数，由 得到与各共振角频率ω_r对应的谐振系统无阻尼时的谐振角频率k，再由得到与各共振角频率ω_r对应的阻尼因子n，其中能量衰减倍数Q为f(t)的振幅的平方与谐振系统与质量无关的能量E’的比例，即定义稳态时表示外强迫力的f(t)的振幅的平方与谐振系统模型因为共振而累积的与质量无关的能量E’的比例定义为能量衰减倍数Q。

通过不同方法得到与各共振角频率ω_r对应的谐振系统的阻尼因子n和谐振系统无阻尼时的谐振角频率k，可根据实际需求选择相应的谐振系统。

当然，可以理解，在本发明的其他实施方式中，也可以使用其他方式来得到谐振系统的k和n。

此后进入步骤103，根据当前采样点对应时刻的状态参数中的速度v和位移x计算谐振强度，作为频域中与共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度。可以理解，在本发明的各个实施方式中，上述谐振强度可以用与质量无关的能量E’=k²x²+v²来表示，也可以用谐振振幅 $A=\sqrt{x^2+\frac{v^2}{k^2}}$ 来表示。

将时域信号用曲线平滑后再生为近似模拟信号，传递给进行有阻尼的受迫简谐振动的谐振系统，并计算该谐振系统的谐振强度，以作为频域信号强度，因为是采用特定形式曲线来作为外强迫力，从而可以直接采用解析解进行精确计算，避免了用诸如欧拉法、威尔逊-θ法或线性加速法等积分算法来求近似解，在提高计算准确性的同时极大地降低了计算量，并且由于可以任选频点计算，所以在只要计算特定的一个或几个频点的信号时，计算量远低于像快速傅里叶变换这类传统计算方法。

更具体地说，该将时域信号转换成频域信号的过程主要如下：

众所周知，对于具有阻尼和外强迫力的弹簧简谐振动系统，定义m为谐振子质量，x为谐振子的瞬时偏移平衡点的位置，v为瞬时速度，a为瞬时加速度。

假设：

1)阻尼作用与瞬时速度v成正比，阻尼系数定为μ，阻力R=-μv；

2)弹簧的弹性系数为K，回复力F=-Kx；

3)外强迫力为Hsinβt，其中H为外力的振幅，β为外力的角频率；

则系统有以下瞬时力平衡关系：ma+μν+Kx＝Hsinβt

将上式以时间微分形式表示： $m\frac{d^{2}x}{{\mathrm{dt}}^2}+\mathrm{\μ}\frac{\mathrm{dx}}{\mathrm{dt}}+\mathrm{Kx}=H\sin \mathrm{\βt}$

为阻尼率；

是谐振系统无阻尼时的固有谐振频率；

设： 其中k就是ω₀，代入上式，则上式可转化为：

$\frac{d^{2}x}{{\mathrm{dt}}^2}+2n\frac{\mathrm{dx}}{\mathrm{dt}}+k^{2}x=h\sin \mathrm{\βt}---\left(1\right)$

这就是有阻尼强迫振动系统(即有阻尼的受迫谐振系统)的去量纲化微分方程。

其中k为谐振系统无阻尼时的固有频率（即谐振角频率），n称为阻尼因子，n/k称为阻尼比，hsinβt称为外强迫力。

可知这个有阻尼谐振系统的无外强迫力时的自由振动角频率ω为：

$\mathrm{\ω}=\sqrt{k^2-n^2}---\left(2\right)$

而这个有阻尼谐振系统的有外强迫力时的共振角频率ω_r为：

${\mathrm{\ω}}_r=\sqrt{k^2-{2n}^2}---\left(3\right)$

众所周知，有阻尼的谐振系统在受到外强迫力f(t)＝hsinβt持续作用后，产生受迫振动的同时，也将外强迫力带来的能量不断的转化为系统自身的动能和势能，同时系统的振幅不断增大。但因为阻尼作用的存在，随着系统振幅的增大，阻尼消耗的能量也逐渐增大，最终系统的传入能量会与阻尼损耗的能量达到平衡，使系统达到一个稳态；而当外强迫力消失，系统的能量会因为阻尼衰减而快速降低，振幅也逐渐减弱到零。

图2展示了这种变化的时间特性。上部曲线为输入信号能量，下部曲线为谐振系统振幅。输入信号的能量反映在谐振系统振幅上，有快速上升，平稳维持和阻尼衰减三个阶段。快速上升阶段T_r，谐振系统接收到的信号能量高于被阻尼衰减掉的能量，因此谐振幅度快速上升。随着谐振振幅逐渐增高，阻尼作用也同时增强，当阻尼所衰减的能量与输入信号的能量达到平衡时，谐振幅度不再增高而是达到一个稳定值。阻尼衰减阶段T_d，当输入信号消失时，随着谐振系统积蓄的能量不断被衰减，谐振振幅也逐渐降低直至消失。

外强迫力的角频率β越接近于谐振系统的共振角频率ω_r，系统的能量积蓄作用就越强，振幅越大。如果两个角频率相等，则谐振幅度将达到最大(即共振状态)。图3展示了这种变化的频率特性，如图3所示，当外强迫力的角频率β等于共振角频率ω_r时，该谐振系统产生最大谐振振幅A_m，而当外强迫力的角频率β偏离共振角频率ω_r，如等于ω_a、ω_b时，该谐振系统的振幅变小为

数字信号（即时域信号）在以时间为横坐标轴，信号幅度为纵坐标轴的波形图上，可以表示为一系列水平间隔相等但垂直位置不同的采样值点，如图4中实心三角形标识，相邻二个采样值点的水平间距就是采样周期；为了从这些采样值点重建原始信号，可以将波形图上每一对相邻采样值点之间用平滑曲线段相连，这些平滑曲线段首尾相接，就构成了近似的连续模拟信号，如图4中连接各采样值点的虚线和实线，每一段平滑曲线可以用f(t)表示，其参数由相邻两个或多个采样值以及必要的其它因素，通过一个与曲线形式对应的算法计算得到。

将(1)式中的外强迫力hsinβt用由离散数据拟合所得的分段平滑曲线f(t)代替，则新的微分方程就可以描述在受到f(t)影响下，谐振系统在一个采样周期中的运动过程，即该系统就构成了一个谐振传感器。

用平滑曲线f(t)代替外强迫力hsinβt后（1）式变为：

$\frac{d^{2}x}{{\mathrm{dt}}^2}+2n\frac{\mathrm{dx}}{\mathrm{dt}}+k^{2}x=f\left(t\right)---\left(4\right)$

对于（4）式，它属于二阶常系数线性非齐次微分方程，它的特征方程 $\frac{d^{2}x}{{\mathrm{dt}}^2}+2n\frac{\mathrm{dx}}{\mathrm{dt}}+k^{2}x=0$ 的根为

因为有阻尼谐振系统的自由振动角频率：

$\mathrm{\ω}={\mathrm{\ω}}_0\sqrt{1-{\mathrm{\ζ}}^2}=\sqrt{\frac{K}{m}-\frac{{\mathrm{\μ}}^2}{{4m}^2}}$

将 代入，得到（2）式，

因为我们要求的是某个实频率ω下的解，这样根据(2)式就有k>n的限制，即阻尼比因此特征方程的根为一对共轭复根：

$r_1=-n+i\sqrt{k^2-n^2},$ $r_2=-n-i\sqrt{k^2-n^2},$ 或表示为：r_1,2＝-n±iω。

所以（4）式的暂态解为：

X＝e^-nt(C₁cosωt+C₂sinωt)   (5)

其中C₁,C₂为待定系数。根据微分方程的求解知识，只要(4)式等号右边的曲线函数f(t)的形式能使微分方程可以得到稳态解y^*，则就可以最终得到x的解析解:

x＝X+y^*   (6)

从而不需要用诸如欧拉法、威尔逊-θ法或线性加速法等积分算法来求近似解，增大计算复杂度。

因为(4)式是二阶常系数非齐次线性方程，对应的可以使用待定系数法，常数变易法以及拉普拉斯变换等方法求解该方程得到解析解。

以待定系数法为例：众所周知，当f(t)满足以下两种常见形式时，我们可以直接得到稳态解y^*：

（1）f(x)＝P_m(x)e^λx，其中λ是常数，P_m(x)是x的一个m次多项式：

P_m(x)＝a₀x^m+a₁x^m-1+…+a_m-1x+a_m；

（2）f(x)＝e^λx[P_l(x)cosωx+P_n(x)sinωx]，其中λ、ω是常数，P_l(x)、P_n(x)分别是x的l次、n次多项式，其中有一个可以为零；

符合上述要求的平滑曲线有多种形式，比如最简单的a+bx直线型，三阶多项式曲线a₀+a₁x+a₂x²+a₃x³型，以及a+be^λx型，b cos px型等。针对上述每一种类型，都是可以得到一个解析解的。

得到了解析解x，也就得到了瞬时位置x和瞬时速度v(即x的导数)的计算公式，就可以由任意一个采样值点所对应时刻的谐振系统的初始位置x₀和初始速度v₀，以及连接下一采样值点的平滑曲线的参数，计算出下一采样值点所对应时刻的瞬时位置x和瞬时速度v。对数字信号连续做上述计算，就可以模拟出谐振系统的振动全过程。

谐振系统的谐振强度可以从每一步计算得到的瞬时位置x和瞬时速度v计算出来。

因为谐振系统积蓄的总能量E等于弹性势能与瞬时动能的和：

$E=\frac{1}{2}{\mathrm{Kx}}^2+\frac{1}{2}{\mathrm{mv}}^2---\left(7\right)$

将代入(7)式，并令则上式转化为与质量m无关的能量值E′＝k²x²+ν²   (8)

因此可以用去质量化的与质量无关的能量E’表示在某一时刻谐振系统的瞬时能量，称为与质量无关的能量。

因为谐振系统的最大位移（即谐振振幅A）出现在瞬时速度v=0的时刻，此时的系统能量和一起代入(7)式后得到谐振振幅A的算法：

$A=\sqrt{x^2+\frac{v^2}{k^2}}---\left(9\right)$

与质量无关的能量E’和谐振振幅A之间的关系为：E’＝k²A²

为了让谐振模型满足不同频率的转换目的，需要根据转换目标来设计谐振模型每一个共振角频率ω_r的参数n和k，这里以hsinβt作为外强迫力来说明参数设计方案。

根据谐振系统的特性，当外强迫力hsinβt的角频率β与谐振系统的共振角频率ω_r相等时，即： ${\mathrm{\ω}}_r=\mathrm{\β}=\sqrt{k^2-{2n}^2}---\left(10\right)$

此时系统有最大谐振振幅A：

第一种谐振模型设计方案：所有共振角频率下的谐振振幅衰减速率相同，即阻尼因子n固定不变；

该设计的特性是谐振振幅的衰减倍数随着固有频率的升高而升高。

设给定的共振角频率为ω_r则按如下步骤确定n和k：

a)给n指定>0的任意实数；

b)由 $k^2={\mathrm{\ω}}_r^2+{2n}^2$ 计算出k。

第二种谐振模型设计方案：要求外强迫力振幅的平方与谐振系统模型因为共振而累积的与质量无关的能量E’的比值对所有频率保持不变；

定义外强迫力h sin βt的振幅的平方h²与与质量无关的能量E’的比值为能量衰减倍数Q，即与(8)（9）式一同代入(11)式,则有如下关系：

4n²(k²-n²)＝k²Q   (12)

因此由(10)式和(12)式可以得到以下设定方法：

设给定的共振角频率为ω_r按如下步骤确定n和k：

a)给Q指定>0的任意实数；

b)由 $k^2=\frac{Q+\sqrt{{4\mathrm{\ω}}_r^4+Q^2}}{2}$ 计算出k²；

c)由 $n=\sqrt{\frac{k^2-{\mathrm{\ω}}_r^2}{2}}$ 计算出n。

当共振角频率ω_r逐渐增大时，k²趋近于n趋近于

该设计的特性是谐振振幅的衰减速度随着频率升高而稳定直至趋近于

第三种谐振模型设计方案：外强迫力振幅与谐振振幅的比例对所有频率保持不变；该设计的特性是谐振振幅的衰减速度随着频率升高而减慢，直至趋近于零。

定义外强迫力h sin βt的振幅与谐振系统模型因为共振而产生的谐振振幅A的比例定义为振幅衰减倍数C，即代入(11)式，则有如下关系：

4n²(k²-n²)＝C²   (13)

因此由(10)式和(13)式可以得到以下设定方法：

设给定的共振角频率为ω_r则按如下步骤确定n和k：

a)给C指定>0的任意实数；

b)由 $k^2=\sqrt{{\mathrm{\ω}}_r^2+C^2}$ 计算出k；

c)由 $n=\sqrt{\frac{k^2-{\mathrm{\ω}}_r^2}{2}}$ 计算出n。

由于瞬时位置x和瞬时速度v是以相邻两个采样值点为区间段进行计算得到的，影响其每一段区间计算结果的因素包括区间的初始状态值x₀,v₀和连接相邻采样值点的表示外强迫力的f(t)的曲线参数，以及与f(t)形式对应的x和v的解析解。所以，在某些特殊需求的场合下，对不同共振频率的谐振模型可以应用不同形式的曲线进行平滑连接；或者更进一步，对同一个共振频率的谐振模型的不同采样值点之间应用不同的表示外强迫力的f(t)的曲线形式，即使用多类型平滑曲线的组合来平滑连接所有采样值点，只要确保始终应用与当前连接相邻采样值点的平滑曲线f(t)形式相对应的x和v的解析解，就可以正确计算x和v，从而实现振动能量累积。

本发明第二实施方式涉及一种将时域信号转换成频域信号的方法。图5是该将时域信号转换成频域信号的方法的流程示意图。

具体地说，如图5所示，第二实施方式在第一实施方式的基础上进行了改进，主要改进之处在于，上述方法还包括以下步骤：

在步骤304中，按时序获取与共振角频率ω_r对应的频点的在时域信号中所有采样点对应时刻的频域信号强度，以形成与共振角频率ω_r对应的频点的时间-频域信号强度曲线。

此后进入步骤305，将预定的频率序列ω₁，…，ω_m乘以2π并依次作为谐振系统的共振角频率ω_r，以得到与各共振角频率ω_r对应的频点的时间-频域信号强度曲线。

此后进入步骤306，将与各共振角频率ω_r对应的频点的时间-频域信号强度曲线按预定的频率序列排列，即得到时间-频率分布图。

通过计算频域信号强度变化过程，进而得到时间-频率分布图，若将相应频点的所有采样点对应时刻的频域信号强度叠加即可得到该频点的频域信号强度，从而成功将时域信号转换成频域信号，因为可以任选频点计算，频率精度与采样频率无关，所以即使时域信号采样频率低，频率精度也可以远高于快速傅里叶变换，而且还可以任意选择感兴趣的频率区间来设置任意多个频率点进行转换，由于是逐个采样点计算，时间响应最高可达到采样周期并且响应幅度可调，因此可以得到同时具有高频率精度和高时间精度的时间-频率分布图。

经过上述改进，转换方法的基本步骤为：

选取一种平滑曲线的函数形式，作为所述谐振系统模型的表示外强迫力的f(t)部分，求解所述谐振系统模型得出所述的瞬时位移x和所述的瞬时速度v的解析解；构造一组能让所述谐振系统模型在指定频率（相当于对应共振角频率ω_r的频点）产生共振的k和n，然后计算出连接数字信号每一对相邻采样值点的平滑曲线段的参数，代入所述的解析解，连续计算所述的指定频率下的每一个所述采样值点对应时刻的所述谐振强度，得到所述指定频率下的时间-谐振强度曲线（即时间-频域信号强度曲线）；计算同一组数字信号在不同所述指定频率下的所述时间-谐振强度曲线，按频率高低顺序排列这些曲线，就得到以谐振强度表示的时间-频率分布图。

作为本发明的一个优选例，以余弦函数作为平滑曲线来进行说明，但不以此为限。

为了使曲线能够在每一个采样值点平滑，并且减少平滑后的波形产生的额外高次谐波，选取余弦函数作为平滑曲线的基本形式，具体形式定为f(t)=a+bcos(pt)，p为数字信号的采样角频率的一半，a为两相邻采样值的算数平均值，b为前一采样值与当前采样值的差的一半。这样相邻的二采样值点就可以用1/2周余弦曲线相连，并且在连接点处的一阶导数连续。如图4中实线所示，三角符号代表采样值，实线代表平滑曲线。

具体计算步骤如下：

i.以单一的a+bcos(pt)形式作为外强迫力，构造一个具有阻尼和强迫力的谐振系统模型，该谐振系统的微分方程为：

$\frac{d^{2}x}{{\mathrm{dt}}^2}+2n\frac{\mathrm{dx}}{\mathrm{dt}}+k^{2}x=a+b\cos \mathrm{pt}---\left(14\right)$

如果我们能够求解上式，就可以通过给出任意一个离散点所对应时刻的谐振系统的初始位移x₀和初始速度v₀，以及连接下一离散点拟合曲线的a,b值，计算出下一离散点所对应时刻的位移x和速度v（即瞬时速度）。

对14式微分方程求解，最终就得到瞬时位置x的解析解，即x的计算公式，和x的一阶导数即瞬时速度v的计算公式

ii.指定一个上述谐振系统的能量衰减倍数Q，Q由实验确定；

iii.设计一个时间-频率分布图的频率序列ω₁，…，ω_m，选取该频率序列的第一个频率点ω₁，并乘以2π作为谐振系统的共振角频率ω_r；

iv.由ω_r计算出k²和n，其中 $k^2=\frac{Q+\sqrt{{4\mathrm{\ω}}_r^4+Q^2}}{2},$ $n=\sqrt{\frac{k^2-{\mathrm{\ω}}_r^2}{2}};$

v.设起始时刻的x₀=0,v₀=0,计算出平滑连接时间零点和数字信号序列中第一个采样值点的平滑曲线段的a,b值，然后用求解(14)式得到的x和v的计算公式计算出第一个采样值点所对应时刻的瞬时位置x和瞬时速度v，再由瞬时位置x和瞬时速度v代入(8)式计算出与质量无关的能量E’或者(9)式计算谐振振幅A作为谐振强度值；

vi.选取数字信号序列的下一个采样值点作为当前采样值点，计算出平滑连接上一个采样值点和当前采样值点的平滑曲线段的a,b值，然后用求解(14)式得到x和v的计算公式计算出当前采样值点所对应时刻的瞬时位置x和瞬时速度v，其中x₀,v₀分别等于前一次计算得到的x和v，再通过(8)式计算出与质量无关的能量E’或者通过(9)式计算谐振振幅A作为谐振强度值；

vii.重复步骤ⅵ直到所有数字信号处理完毕，将得到的所有谐振强度绘制在以时间为x轴，谐振强度为y轴的二维平面上，得到数字信号在指定谐振频率下的时间-谐振强度曲线；

viii.从步骤ⅲ的频率序列中选取下一个频率点ω_i，并乘以2π作为谐振系统的共振角频率ω_r，执行步骤ⅳ～ⅶ，计算出新的频率点ω_i下的时间-谐振强度曲线；

ix.重复执行步骤ⅷ直到ⅲ的频率序列中所有频率点的时间-谐振强度曲线都计算完毕，将所有曲线按频率顺序排列，就得到了以谐振强度表示的时间-频率分布图。

通过连续计算所有的离散点，这样我们就得到了在某一个谐振频率点上的能量幅度-时间响应曲线。

对任意选取的一系列频率点构建相应谐振频率的谐振系统，计算同一组离散信号，并将所产生的能量幅度-时间响应曲线按频率进行合并，就构成了三维谐振时间-频率分布图。

本发明第三实施方式涉及一种将时域信号转换成频域信号的方法。图6是该将时域信号转换成频域信号的方法的流程示意图。

具体地说，如图6所示，第三实施方式在第一实施方式的基础上进行了改进，主要改进之处在于，上述方法还包括以下步骤：

在步骤404中，将预定的频率序列ω₁，…，ω_m乘以2π并依次作为谐振系统的共振角频率ω_r，得到与各共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度，以形成在当前采样点对应时刻的频率-频域信号强度曲线。

此后进入步骤405，按时序获取时域信号中所有采样点对应时刻的频率-频域信号强度曲线。

此后进入步骤406，将所有采样点对应时刻的频率-频域信号强度曲线按时序排列，即得到时间-频率分布图。

通过计算频域信号强度变化过程，进而得到时间-频率分布图，若将相应频点的所有采样点对应时刻的频域信号强度叠加即可得到该频点的频域信号强度，从而成功将时域信号转换成频域信号，因为可以任选频点计算，频率精度与采样频率无关，所以即使时域信号采样频率低，频率精度也可以远高于快速傅里叶变换，而且还可以任意选择感兴趣的频率区间来设置任意多个频率点进行转换，由于是逐个采样点计算，时间响应最高可达到采样周期并且响应幅度可调，因此可以得到同时具有高频率精度和高时间精度的时间-频率分布图。

经过上述改进，转换方法的基本步骤为：

先根据平滑曲线的函数形式推导出所述谐振系统模型相应的解析解，并预先建立时间-频率分布图频率序列的所有频率点相对应的所述谐振系统模型的k和n；对一个采样值点先计算连接前一个采样值点的所述平滑曲线的参数，再用这个参数计算所有频率点下的所述谐振强度，然后再对下一个采样值点做同样计算直至所有采样值点计算完毕，最终同时得到所有频率点的所述时间-谐振强度曲线，再汇集成所述以谐振强度表示的时间-频率分布图。

作为本发明的一个优选例，具体计算步骤以单一的a′+b′p′t形式为例，但不以此为限。

在真实信号频率远低于采样频率，并且噪音较低的场合，可以使用简单的直线a+bx，或者更精确的三阶多项式曲线a₀+a₁x+a₂x²+a₃x³来平滑连接采样值点，从而降低额外产生的高次谐波。

对于直线形式f(t)=a′+b′p′t，其参数的确定方法是：p′是数字信号的采样频率，a′取前一采样值，b′取前一采样值与当前采样值的差值。

对于三阶多项式形式f(t)＝a₀+a₁p′t+a₂p′²t²+a₃p′³t³，其参数的确定方法可以由3次贝赛尔曲线控制点的确定方法推导而得到，其中p′是数字信号的采样频率。

为了减少计算量，可以对每一对采样值点，先计算平滑曲线段的参数，然后用这段曲线计算所有共振频率模型的瞬时位置x和瞬时速度v。

i.以a′+b′p′t形式作为外强迫力，构造一个具有阻尼和强迫力的简谐振动系统模型，该简谐振动系统的微分方程为：

$\frac{d^{2}x}{{\mathrm{dt}}^2}+2n\frac{\mathrm{dx}}{\mathrm{dt}}+k^{2}x=a^{\′}+-b^{\′}{p}^{\′}t--\left(17\right)$

求解过程与第三实施方式中的优选例的步骤ⅰ相同；

ii.指定一个上述谐振系统的阻尼因子n，n由实验确定；

iii.设计一个时间-频率分布图的频率序列ω₁，…，ω_m，由计算出所有频率点所对应的k²；

iv.指定t=0时的x₀=0,v₀=0作为初始的前一采样值点；

v.从数字信号序列中按时间顺序取一个新采样值点作为当前采样值点，计算平滑连接前一采样值点和当前采样值点的平滑曲线段的参数值；

vi.然后用求解(17)式得到的x和v的计算公式，分别对频率序列中每一个频率点，计算出当前采样值点所对应时刻的瞬时位置x和瞬时速度v，再由瞬时位置x和瞬时速度v代入(8)式计算出与质量无关的能量E’或者(9)式计算谐振振幅A作为谐振强度值；

vii.将得到的所有频率点的谐振强度绘制在以频率为x轴谐振强度为y轴的二维平面上，得到数字信号在第一个采样值点所对应时刻的频率-谐振强度曲线（即频率-频域信号强度曲线）；

viii.重复执行步骤ⅴ～ⅶ，直到所有数字信号采样值都计算完毕，将所有频率-谐振强度曲线按时间顺序排列，就得到了以谐振强度表示的时间-频率分布图。

本发明第四实施方式涉及一种声纹识别方法。该声纹识别方法在声纹识别的过程中使用第一、第二或第三实施方式中的方法根据声音信号产生时间-频率分布图。

本实施方式能提高声纹识别的准确度，因为本方法对声音信号所产生的时频图可以比别的方法得到的更为精细。

作为本发明的优选例，选取一些有代表性的声音片段，如语言片段、粉红噪声、调频调幅混合声等进行转换，由短时傅里叶变换和本发明转换得到的频谱图结果如图7a至7g所示。其中：

图7a为“goodbye”语音片段的STFT(短时傅里叶频谱图)，图像在时间轴上放大4倍，在频率轴上放大4倍。Sample Rate:11025Hz，SampleTime:0.76s，Window Type:Hann Window，Window Size:128，FFT FrameSize:256。帧重叠率:50%

图7b和图7c为“goodbye”语音片段的谐振频谱图（即所述时间-频率分布图）,图像在时间轴上压缩20倍，Sample Rate=11025Hz,SampleTime=0.76s,频率点分布类型为对数分布，频率分割精度(点)=512,最低显示频率=64Hz。在图7b中，阻尼因子n=100；而在图7c中，阻尼因子n=50。图7b和图7c显示了能量衰减时产生的拖尾峰。

图7d为粉红噪声片段的STFT(短时傅里叶频谱图),图像在时间轴上放大4倍，在频率轴上放大4倍。Sample Rate:44100Hz，Sample Time:0.14s，Window Type:Hann Window，Window Size:128，FFT Frame Size:256。

图7e为粉红噪声片段的谐振频谱图,图像在时间轴压缩80倍。SampleRate=44100Hz,Sample Time=0.14s,频率分割精度(点)=512。频率点分布类型=线性分布,最低显示频率=64Hz,阻尼因子n=100。

图7f为调频调幅片段的STFT(时频图),图像在时间轴上放大4倍，在频率轴上放大4倍。Sample Rate:5000Hz，Sample Time:0.81s，Window Type:Hann Window，Window Size:128，FFT Frame Size:256。

图7g为调频调幅片段的谐振频谱图,图像在时间轴压缩20倍。SampleRate=5000Hz,Sample Time=0.81s,频率分割精度(点)=512。频率点分布类型＝线性分布，最低显示频率＝50Hz，阻尼因子n＝50。

试验结果归一化为256级灰度图像，谐振频谱图（即所述时间-频率分布图）图像的每一点取时间轴20个点的平均值。

由图7a～图7g可以得出下面的结论：

（1）一个最明显的特征是所有的峰都有拖尾，这与谐振系统能量积蓄和衰减的过程是一致的。

（2）当试验设计的谐振系统的阻尼因子变小时(即具有大的品质因数)，拖尾将变得更长，但波峰在频率方向的分布变得更窄，使频率分辨率提高。

（3）对于真实的语音信号，所产生的谐振频谱图具有比用STFT产生的频谱图高得多的清晰度，不论是在时间响应上还是在频率分辨率上。

本发明第五实施方式涉及一种地震波分析方法。该地震波分析方法使用第一、第二或第三实施方式中的方法根据地震波产生时间-频率分布图，以对地震波进行时频分析。

由于本方法可以对低频率的地震波进行高精度和高时变性的分析，能得出更为精确的地质断层分析结果。

实时对地震波进行高精度的时频分析，有助于提高对地震发生前地质断层活动模式的辨识精度，从而提高地震预测水平。

本发明第六实施方式涉及一种电力线路谐波分析方法。该电力线路谐波分析方法使用第一、第二或第三实施方式中的方法根据电力线路谐波产生时间-频率分布图，以对电力线路谐波进行分析。

由于没有FFT产生的频谱泄露缺陷以及栅栏效应，分析的准确度和时间响应度都会比现存方法高。

本发明第七实施方式涉及一种雷达的回波分析方法。该雷达的回波分析方法使用第一、第二或第三实施方式中的方法根据雷达的回波信号产生时间-频率分布图，以对回波信号进行时频分析。

声纳和雷达的回波一般都是具有时变特性的，传统的傅里叶变换方法以及小波变换方法仍然难以对时变信号做频率精度和时间精度同时都要求很高的分析。本发明的方法可以同时在时间和频率两个维度更精确的分析回波信号，从而提高对探测目标特征的辨识精度。

该时域转换成频域的方法除了可以运用在上述地质勘探，声纹识别，地震波分析，电力线路谐波分析，雷达的回波分析的领域中，还可以运用于：

语音识别，由于频率精度可以划分的很细，可以直接得到某一个时间点所有共振峰的频率点，由此可以精确地推算出基频，并得出准确的共振峰包络线，有助于提高语音识别的准确度；

机械振动分析，用传统的方法分析机械运转产生的振动，只能对振动的平均频率分布做分析，难以对随时间快速变化的振动声音（如低速柴油机的塔塔声）进行精确的分析，而利用本方法对时变变信号的分析能力，可以使之变为可能，从而提高对机械运转故障前兆的预测能力；

人体心脏，大脑，骨骼运动，血液流动等领域的信号分析；

超声波诊断，利用的是超声波声纳技术。等等。

本发明第八实施方式涉及一种将时域信号转换成频域信号的系统。图8是该将时域信号转换成频域信号的系统的结构示意图。

具体地说，如图8所示，该将时域信号转换成频域信号的系统包括：

平滑模块，用于根据时域信号中前一采样点和当前采样点确定平滑曲线f(t)的函数形式和相应参数。

谐振模块，用于将f(t)的参数和前一采样点对应时刻的状态参数输入到一个偏微分方程的解析解中，求得当前采样点对应时刻的状态参数，其中该偏微分方程代表一个有阻尼的受迫谐振系统，形式为 k是该谐振系统无阻尼时的谐振角频率，n为该谐振系统的阻尼因子，k和n对应该谐振系统的一个共振角频率ω_r，状态参数中包含该谐振系统的速度v和位移x。

计算模块，用于根据当前采样点对应时刻的状态参数中的速度v和位移x计算谐振强度，作为频域中与共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度。可以理解，在本发明的各个实施方式中，上述谐振强度可以用与质量无关的能量E’=k²x²+v²来表示，也可以用谐振振幅 $A=\sqrt{x^2+\frac{v^2}{k^2}}$ 来表示。

该系统包括平滑模块、谐振模块和计算模块，其中该平滑模块将时域信号用曲线平滑后再生生近似模拟信号，传递给该谐振模块进行有阻尼的受迫简谐振动，计算模块计算该受迫简谐振动的谐振强度，以作为频域信号强度，由于可以任选频点计算，所以在只要计算特定的一个或几个频点的信号时，计算量远低于传统计算方法，并且由于是可以直接采用解析解进行计算，从而不需要用诸如欧拉法、威尔逊-θ法或线性加速法等积分算法来求近似解，在保持计算准确性的同时极大地降低了计算复杂度。

可以理解，平滑模块可以用于对于时域信号中各采样点和具有不同共振角频率ω_r的受谐振系统，将f(t)设为相同平滑曲线函数形式，也可以用于对于时域信号中各采样点或具有不同共振角频率ω_r的谐振系统，将f(t)设为不同平滑曲线函数形式。

例如，平滑模块可以用于将f(t)设为a+b cos(pt)的形式，其中p为时域信号的采样角频率的一半，a为前一采样值与当前采样值的算数平均值，b为前一采样值与当前采样值的差的一半，可以用于将f(t)设为a′+b′p′t的形式，其中p′为时域信号的采样频率，a′为前一采样值，b′为前一采样值与当前采样值的差，可以用于将f(t)设为a₀+a₁p′t+a₂p′²t²+a₃p′³t³的形式，其中p′为时域信号的采样频率，也可以用于将f(t)设为a+bcos(pt)、a’+b’p’t和a₀+a₁p′t+a₂p’²t²+a₃p’³t³的形式组合。

此外，可以理解，在本发明的其他实施方式中，f(t)还可以采用其他的平滑曲线函数形式。

平滑模块对时域信号中各采样点或具有不同共振角频率的谐振系统，将f(t)设为相同或不同的曲线函数形式进行平滑连接，可适应一些特殊需求的场合。

此外，在本发明的各个实施方式中，谐振模块可以通过以下几个方式来得到谐振系统的k和n：

1）谐振模块固定谐振系统的阻尼因子n不变且设n为大于零的实数，并由得到与各共振角频率ω_r对应的谐振系统无阻尼时的谐振角频率k；

2）谐振模块固定谐振系统的振幅衰减倍数C不变且设C为大于零的实数，并由得到与各共振角频率ω_r对应的谐振系统无阻尼时的谐振角频率k，再由得到与各共振角频率ω_r对应的阻尼因子n，其中振幅衰减倍数C为f(t)的振幅与谐振系统的谐振振幅A的比例。

3）谐振模块固定谐振系统的能量衰减倍数Q不变且设Q为大于零的实数，并由得到与各共振角频率ω_r对应的谐振系统无阻尼时的谐振角频率k，再由得到与各共振角频率ω_r对应的阻尼因子n，其中能量衰减倍数Q为f(t)的振幅的平方与谐振系统与质量无关的能量E′的比例。

谐振模块通过不同方法得到与各共振角频率ω_r对应的谐振系统的阻尼因子n和谐振系统无阻尼时的谐振角频率k，可根据实际需求选择相应的谐振系统。

当然，可以理解，在本发明的其他实施方式中，也可以使用其他方式来得到谐振系统的k和n。

第一实施方式是与本实施方式相对应的方法实施方式，本实施方式可与第一实施方式互相配合实施。第一实施方式中提到的相关技术细节在本实施方式中依然有效，为了减少重复，这里不再赘述。相应地，本实施方式中提到的相关技术细节也可应用在第一实施方式中。

本发明第九实施方式涉及一种将时域信号转换成频域信号的系统。图9是该将时域信号转换成频域信号的系统的结构示意图。

具体地说，如图9所示，第九实施方式在第八实施方式的基础上进行了改进，主要改进之处在于：

该系统还包括构建模块。

上述计算模块用于按时序获取与共振角频率ω_r对应的频点的在时域信号中所有采样点对应时刻的频域信号强度，并按预定的频率序列ω₁，…，ω_m获取与各共振角频率ω_r对应的频点的在时域信号中所有采样点对应时刻的频域信号强度。

谐振模块用于将预定的频率序列ω₁，…，ω_m乘以2π并依次作为谐振系统的共振角频率ω_r。

构建模块用于根据与各共振角频率ω_r对应的频点的在时域信号中所有采样点对应时刻的频域信号强度构建与各共振角频率ω_r对应的频点的时间-频域信号强度曲线，并且将与各共振角频率ω_r对应的频点的时间-频域信号强度曲线按预定的频率序列排列，以得到时间-频率分布图。

通过计算模块计算频域信号强度变化过程，进而由构建模块得到时间-频率分布图，若再将相应频点的所有采样点对应时刻的频域信号强度叠加即可得到该频点的频域信号强度，从而成功将时域信号转换成频域信号，因为可以任选频点计算，频率精度与采样频率无关，所以即使时域信号采样频率低，频率精度也可以远高于快速傅里叶变换，而且还可以任意选择感兴趣的频率区间来设置任意多个频率点进行转换，由于是逐个采样点计算，时间响应最高可达到采样周期并且响应幅度可调，因此可以得到同时具有高频率精度和高时间精度的时间-频率分布图。

第二实施方式是与本实施方式相对应的方法实施方式，本实施方式可与第二实施方式互相配合实施。第二实施方式中提到的相关技术细节在本实施方式中依然有效，为了减少重复，这里不再赘述。相应地，本实施方式中提到的相关技术细节也可应用在第二实施方式中。

本发明第十实施方式涉及一种将时域信号转换成频域信号的系统。

第十实施方式与第九实施方式基本相同，区别主要在于，在第九实施方式中：

计算模块用于按时序获取与共振角频率ω_r对应的频点的在时域信号中所有采样点对应时刻的频域信号强度，并按预定的频率序列ω₁，…，ω_m获取与各共振角频率ω_r对应的频点的在时域信号中所有采样点对应时刻的频域信号强度。

谐振模块用于将预定的频率序列ω₁，…，ω_m乘以2π并依次作为谐振系统的共振角频率ω_r。

构建模块用于根据与各共振角频率ω_r对应的频点的在时域信号中所有采样点对应时刻的频域信号强度构建与各共振角频率ω_r对应的频点的时间-频域信号强度曲线，并且将与各共振角频率ω_r对应的频点的时间-频域信号强度曲线按预定的频率序列排列，以得到时间-频率分布图。

然而在第十实施方式中：

计算模块用于按预定的频率序列ω₁，…，ω_m获取与各共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度，并按时序获取与各共振角频率ω_r对应的频点的在时域信号中所有采样点对应时刻的频域信号强度。

谐振模块用于将预定的频率序列ω₁，…，ω_m乘以2π并依次作为谐振系统的共振角频率ω_r。

构建模块用于根据与各共振角频率ω_r对应的频点的在时域信号中所有采样点对应时刻的频域信号强度构建时域信号中所有采样点对应时刻的频率-频域信号强度曲线，并且将所有采样点对应时刻的频率-频域信号强度曲线按时序排列，以得到时间-频率分布图。

通过计算模块计算频域信号强度变化过程，进而由构建模块得到时间-频率分布图，若再将相应频点的所有采样点对应时刻的频域信号强度叠加即可得到该频点的频域信号强度，从而成功将时域信号转换成频域信号，因为可以任选频点计算，频率精度与采样频率无关，所以即使时域信号采样频率低，频率精度也可以远高于快速傅里叶变换，而且还可以任意选择感兴趣的频率区间来设置任意多个频率点进行转换，由于是逐个采样点计算，时间响应最高可接近采样周期，并且响应幅度可调，因此可以得到同时具有高频率精度和高时间精度的时间-频率分布图。

第三实施方式是与本实施方式相对应的方法实施方式，本实施方式可与第三实施方式互相配合实施。第三实施方式中提到的相关技术细节在本实施方式中依然有效，为了减少重复，这里不再赘述。相应地，本实施方式中提到的相关技术细节也可应用在第三实施方式中。

本发明采用模拟有阻尼的简谐振动系统振动过程的方法，来计算数字信号在某一特定频率下的谐振强度变化过程。具体做法是：通过一系列平滑曲线段平滑连接数字信号在波形图上的每一个采样值点，从而将数字信号转换为近似的连续模拟信号；将这些平滑曲线段作为外强迫力依次代入一个具有某个共振频率的模拟谐振系统模型进行有阻尼的受迫振动，通过计算每一段平滑曲线段所导致的谐振强度变化，最终得到谐振强度的变化过程即时间-谐振强度曲线；同样的，让同一组平滑曲线代入一组各自具有不同共振频率的模拟谐振系统模型进行计算，就可以得到一组不同频率下的时间-谐振强度曲线，再按频率顺序排列这些曲线，就生成了以谐振强度表示的时间-频率分布图。

本发明针对短时傅里叶变换计算时间频率分布图结果精度低的缺陷，提出了一种将数字信号用曲线平滑再生为近似模拟信号，传递给模拟谐振系统模型做有阻尼的受迫简谐振动，通过计算谐振强度变化过程，进而计算出时间频率分布图的算法。该时间频率分布图与短时傅里叶变换的结果在形式上类似，并也具有真实的物理意义。该算法的优点是时间响应可调，并且精度高，最高可以接近采样周期；同时，共振频率可以在低于奈奎斯特频率的范围内任意指定，频率准确度和精度高，并且时间精度和频率精度之间没有互相制约关系。

需要说明的是，本发明各设备实施方式中提到的各单元都是逻辑单元，在物理上，一个逻辑单元可以是一个物理单元，也可以是一个物理单元的一部分，还可以以多个物理单元的组合实现，这些逻辑单元本身的物理实现方式并不是最重要的，这些逻辑单元所实现的功能的组合才是解决本发明所提出的技术问题的关键。此外，为了突出本发明的创新部分，本发明上述各设备实施方式并没有将与解决本发明所提出的技术问题关系不太密切的单元引入，这并不表明上述设备实施方式并不存在其它的单元。

需要说明的是，在本专利的权利要求和说明书中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

虽然通过参照本发明的某些优选实施例，已经对本发明进行了图示和描述，但本领域的普通技术人员应该明白，可以在形式上和细节上对其作各种改变，而不偏离本发明的精神和范围。

Claims (18)

1.一种将时域信号转换成频域信号的方法，其特征在于，包括以下步骤：

根据时域信号中前一采样点和当前采样点确定平滑曲线f(t)的函数形式和相应参数；

将所述f(t)的参数和前一采样点对应时刻的状态参数输入到一个偏微分方程的解析解中，求得当前采样点对应时刻的状态参数，其中该偏微分方程代表一个有阻尼的受迫谐振系统，形式为f(t)是所述平滑曲线，k是该谐振系统无阻尼时的谐振角频率，n为该谐振系统的阻尼因子，k和n对应该谐振系统的一个共振角频率ω_r，状态参数中包含该谐振系统的速度v和位移x；

根据当前采样点对应时刻的状态参数中的速度v和位移x计算谐振强度，作为频域中与所述共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度。

2.根据权利要求1所述的将时域信号转换成频域信号的方法，其特征在于，还包括以下步骤：

按时序获取与所述共振角频率ωr对应的频点的在所述时域信号中所有采样点对应时刻的频域信号强度，以形成与所述共振角频率ω_r对应的频点的时间－频域信号强度曲线；

将预定的频率序列ω₁，…，ω_m乘以2π并依次作为所述谐振系统的共振角频率ω_r，以得到与各共振角频率ω_r对应的频点的时间－频域信号强度曲线；

将与各共振角频率ω_r对应的频点的所述时间－频域信号强度曲线按预定的频率序列排列，即得到时间－频率分布图。

3.根据权利要求1所述的将时域信号转换成频域信号的方法，其特征在于，还包括以下步骤：

将预定的频率序列ω₁，…，ω_m乘以2π并依次作为所述谐振系统的共振角频率ω_r，得到与各共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度，以形成在当前采样点对应时刻的频率－频域信号强度曲线；

按时序获取所述时域信号中所有采样点对应时刻的频率－频域信号强度曲线；

将所有采样点对应时刻的频率－频域信号强度曲线按时序排列，即得到时间－频率分布图。

4.根据权利要求1至3中任一项所述的将时域信号转换成频域信号的方法，其特征在于，所述f(t)对于所述时域信号中各采样点和具有不同共振角频率ω_r的谐振系统采用相同平滑曲线函数形式；

或所述f(t)对于所述时域信号中各采样点采用不同平滑曲线函数形式；

或所述f(t)对于具有不同共振角频率ω_r的谐振系统采用不同平滑曲线函数形式。

5.根据权利要求4所述的将时域信号转换成频域信号的方法，其特征在于，所述f(t)采用a+bcos(pt)的形式，其中p为所述时域信号的采样角频率的一半，a为前一采样值与当前采样值的算数平均值，b为前一采样值与当前采样值的差的一半；

或所述f(t)采用a′+b′p′t的形式，其中p′为所述时域信号的采样频率，a′为前一采样值，b′为前一采样值与当前采样值的差；

或所述f(t)采用a₀+a₁p′t+a₂p′²t²+a₃p′³t³的形式，其中p′为所述时域信号的采样频率；

或所述f(t)采用a+bcos(pt)、a′+b′p′t和a₀+a₁p′t+a₂p′²t²+a₃p′³t³的形式组合。

6.根据权利要求1至3中任一项所述的将时域信号转换成频域信号的方法，其特征在于，所述谐振强度用与质量无关的能量E′=k²x²+v²来表示；或所述谐振强度用谐振振幅来表示。

7.根据权利要求6所述的将时域信号转换成频域信号的方法，其特征在于，所述谐振系统的阻尼因子n不变且为大于零的实数，由得到与各共振角频率ω_r对应的所述谐振系统无阻尼时的谐振角频率k；

或所述谐振系统的振幅衰减倍数C不变且为大于零的实数，由得到与各共振角频率ω_r对应的所述谐振系统无阻尼时的谐振角频率k，再由得到与各共振角频率ω_r对应的阻尼因子n，其中振幅衰减倍数C为所述f(t)的振幅与所述谐振系统的谐振振幅A的比例；

或所述谐振系统的能量衰减倍数Q不变且为大于零的实数，由得到与各共振角频率ω_r对应的所述谐振系统无阻尼时的谐振角频率k，再由得到与各共振角频率ω_r对应的阻尼因子n，其中能量衰减倍数Q为所述f(t)的振幅的平方与所述谐振系统与质量无关的能量E′的比例。

8.一种声纹识别方法，其特征在于，在声纹识别的过程中使用权利要求1所述的方法根据声音信号产生时间－频率分布图。

9.一种地震波分析方法，其特征在于，使用权利要求1所述的方法根据地震波产生时间－频率分布图，以对地震波进行时频分析。

10.一种电力线路谐波分析方法，其特征在于，使用权利要求1所述的方法根据电力线路谐波产生时间－频率分布图，以对电力线路谐波进行分析。

11.一种雷达的回波分析方法，其特征在于，使用权利要求1所述的方法根据雷达的回波信号产生时间－频率分布图，以对回波信号进行时频分析。

12.一种将时域信号转换成频域信号的系统，其特征在于，包括：

平滑模块，用于根据时域信号中前一采样点和当前采样点确定平滑曲线f(t)的函数形式和相应参数；

谐振模块，用于将所述f(t)的参数和前一采样点对应时刻的状态参数输入到一个偏微分方程的解析解中，求得当前采样点对应时刻的状态参数，其中该偏微分方程代表一个有阻尼的受迫谐振系统，形式为 k是该谐振系统无阻尼时的谐振角频率，n为该谐振系统的阻尼因子，k和n对应该谐振系统的一个共振角频率ω_r，状态参数中包含该谐振系统的速度v和位移x；

计算模块，用于根据当前采样点对应时刻的状态参数中的速度v和位移x计算谐振强度，作为频域中与所述共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度。

13.根据权利要求12所述的将时域信号转换成频域信号的系统，其特征在于，还包括构建模块；

所述计算模块用于按时序获取与所述共振角频率ω_r对应的频点的在所述时域信号中所有采样点对应时刻的频域信号强度，并按预定的频率序列ω₁，…，ω_m获取与各共振角频率ω_r对应的频点的在所述时域信号中所有采样点对应时刻的频域信号强度；

所述谐振模块用于将预定的频率序列ω₁，…，ω_m乘以2π并依次作为所述谐振系统的共振角频率ω_r；

所述构建模块用于根据与各共振角频率ω_r对应的频点的在所述时域信号中所有采样点对应时刻的频域信号强度构建与各共振角频率ω_r对应的频点的时间－频域信号强度曲线，并且将与各共振角频率ω_r对应的频点的所述时间－频域信号强度曲线按预定的频率序列排列，以得到时间－频率分布图。

14.根据权利要求12所述的将时域信号转换成频域信号的系统，其特征在于，还包括构建模块；

所述计算模块用于按预定的频率序列ω₁，…，ω_m获取与各共振角频率ω_r对应的频点的在当前采样点对应时刻的频域信号强度，并按时序获取与各共振角频率ω_r对应的频点的在所述时域信号中所有采样点对应时刻的频域信号强度；

所述谐振模块用于将预定的频率序列ω₁，…，ω_m乘以2π并依次作为所述谐振系统的共振角频率ω_r；

所述构建模块用于根据与各共振角频率ω_r对应的频点的在所述时域信号中所有采样点对应时刻的频域信号强度构建所述时域信号中所有采样点对应时刻的频率－频域信号强度曲线，并且将所有采样点对应时刻的频率－频域信号强度曲线按时序排列，以得到时间－频率分布图。

15.根据权利要求12至14中任一项所述的将时域信号转换成频域信号的系统，其特征在于，所述平滑模块用于对于所述时域信号中各采样点和具有不同共振角频率ω_r的谐振系统，将所述f(t)设为相同平滑曲线函数形式；

或所述平滑模块用于对于所述时域信号中各采样点，将所述f(t)设为不同平滑曲线函数形式；

或所述平滑模块用于对于具有不同共振角频率ω_r的谐振系统，将所述f(t)设为不同平滑曲线函数形式。

16.根据权利要求15所述的将时域信号转换成频域信号的系统，其特征在于，所述平滑模块用于将所述f(t)设为a+bcos(pt)的形式，其中p为所述时域信号的采样角频率的一半，a为前一采样值与当前采样值的算数平均值，b为前一采样值与当前采样值的差的一半；

或所述平滑模块用于将f(t)设为a′+b′p′t的形式，其中p′为所述时域信号的采样频率，a′为前一采样值，b′为前一采样值与当前采样值的差；

或所述平滑模块用于将所述f(t)设为a₀+a₁p′t+a₂p′²t²+a₃p′³t³的形式，其中p′为所述时域信号的采样频率；

或所述平滑模块用于将所述f(t)设为a+bcos(pt)、a′+b′p′t和a₀+a₁p′t+a₂p′²t²+a₃p′³t³的形式组合。

17.根据权利要求12至14中任一项所述的将时域信号转换成频域信号的系统，其特征在于，所述谐振强度用与质量无关的能量E′=k²x²+v²来表示；

或所述谐振强度用谐振振幅来表示。

18.根据权利要求17所述的将时域信号转换成频域信号的系统，其特征在于，所述谐振模块用于固定所述谐振系统的阻尼因子n不变且设n为大于零的实数，并由得到与各共振角频率ω_r对应的所述谐振系统无阻尼时的谐振角频率k；

或所述谐振模块用于固定所述谐振系统的振幅衰减倍数C不变且设C为大于零的实数，并由得到与各共振角频率ω_r对应的所述谐振系统无阻尼时的谐振角频率k，再由得到与各共振角频率ω_r对应的阻尼因子n，其中振幅衰减倍数C为所述f(t)的振幅与所述谐振系统的谐振振幅A的比例；

或所述谐振模块用于固定所述谐振系统的能量衰减倍数Q不变且设Q为大于零的实数，并由得到与各共振角频率ω_r对应的所述谐振系统无阻尼时的谐振角频率k，再由得到与各共振角频率ω_r对应的阻尼因子n，其中能量衰减倍数Q为所述f(t)的振幅的平方与所述谐振系统与质量无关的能量E′的比例。