CN106415310A - 峰值频率检测装置、方法以及程序 - Google Patents

峰值频率检测装置、方法以及程序 Download PDF

Info

Publication number
CN106415310A
CN106415310A CN201480079481.8A CN201480079481A CN106415310A CN 106415310 A CN106415310 A CN 106415310A CN 201480079481 A CN201480079481 A CN 201480079481A CN 106415310 A CN106415310 A CN 106415310A
Authority
CN
China
Prior art keywords
frequency
digital data
power
data strings
band
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201480079481.8A
Other languages
English (en)
Other versions
CN106415310B (zh
Inventor
阿久津史
阿久津一史
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Kk Sfft
Original Assignee
Kk Sfft
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Kk Sfft filed Critical Kk Sfft
Publication of CN106415310A publication Critical patent/CN106415310A/zh
Application granted granted Critical
Publication of CN106415310B publication Critical patent/CN106415310B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R19/00—Arrangements for measuring currents or voltages or for indicating presence or sign thereof
    • G01R19/04—Measuring peak values or amplitude or envelope of ac or of pulses
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis
    • G01R23/165—Spectrum analysis; Fourier analysis using filters
    • G01R23/167—Spectrum analysis; Fourier analysis using filters with digital filters
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis
    • G01R23/177—Analysis of very low frequencies
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/02—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems using reflection of acoustic waves
    • G01S15/50—Systems of measurement, based on relative movement of the target
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/35—Details of non-pulse systems
    • G01S7/352—Receivers
    • G01S7/356—Receivers involving particularities of FFT processing
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/523—Details of pulse systems
    • G01S7/526—Receivers
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10G—REPRESENTATION OF MUSIC; RECORDING MUSIC IN NOTATION FORM; ACCESSORIES FOR MUSIC OR MUSICAL INSTRUMENTS NOT OTHERWISE PROVIDED FOR, e.g. SUPPORTS
    • G10G7/00—Other auxiliary devices or accessories, e.g. conductors' batons or separate holders for resin or strings
    • G10G7/02—Tuning forks or like devices
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/90—Pitch determination of speech signals

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Acoustics & Sound (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Multimedia (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Human Computer Interaction (AREA)
  • Signal Processing (AREA)
  • Computational Linguistics (AREA)
  • Discrete Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Power Engineering (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Measuring Frequencies, Analyzing Spectra (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明能避免存在于频率分辨率(f0)与时间窗长度(T0)之间的互反关系，能以所希望的频率分辨率和所希望的时间窗长度来检测信号波的峰值频率。本发明提供一种峰值频率检测装置，检测规定的频带(fcl～fch)中功率谱最大的峰值频率，具备：n次方部，对数字数据串的各元素取n次方(n是大于等于2的整数)；FFT部，将与对根据采样频率fs、频率分辨率ftg、时间窗长度Ttg确定的N(N是2的乘方的整数)个取n次方后的采样频率fs的数字数据串进行快速傅立叶转换所得的功率谱的最大值对应的频率作为虚拟峰值频率导出；以及1/n倍部，将对所述虚拟峰值频率取1/n倍后的值作为数字数据串的峰值频率输出，所述峰值频率检测装置满足：n≥1/(ftg×Ttg)、fs/(n×ftg)≤N≤fs×Ttg、fs＞2×n×fch。

Description

峰值频率检测装置、方法以及程序
技术领域
本发明涉及功率谱最大的峰值频率的检测装置、方法以及程序。
背景技术
当将快速傅立叶转换(FFT)的时间窗长度设为T0(s)、将采样频率设为fs(Hz)、将采样数设为N(N是2的乘方的整数)时，则存在如下关系：
[算式1]
T0＝N/fs。
时间窗长度的倒数f0(Hz)称为频率分辨率，为如下关系：
[算式2]
f0＝1/T0。
在通过FFT对信号波进行频率分析的情况下，时间窗长度是进行FFT的信号波的时间长度，频率分辨率是信号波频率的最小分辨率，即为频率的检测精度。以下，也将该时间窗长度称为“进行FFT的信号波的截取时间”、“进行FFT的信号波的时间”等。[算式2]的意思是：FFT的时间窗长度与频率分辨率呈互反关系。该互反关系在频率分辨率小的频率处影响显著。例如，如果需要0.01(Hz)的频率分辨率，则必须对其倒数的100(s)期间的信号波进行FFT。此外，如果进行FFT的信号波的时间是0.01(s)，则只能得到其倒数100(Hz)的频率分辨率。
当利用FFT时，该互反关系有时成为制约，为了避免互反关系的制约，想尽了各种办法。例如，可想到如下方法：在对接收的信号波进行FFT，求出构成在频域出现的频谱的峰值频率时，将时间窗长度T0固定为预定值，对从所接收的信号波以间隔T0而不相互重叠的方式错开并截取的多个期间分别应用FFT(即，多次应用FFT)，将求得的多个峰值频率取平均。在该方法中，有时也能以高于通过一次FFT求得的频率分辨率f0的分辨率(细频率间隔)计算出峰值频率。但是，存在如下问题：如果增加FFT的操作次数，则进行FFT的信号波的截取时间随着其次数增加而变长。
需要说明的是，当对数字数据串应用FFT时，对于提供给FFT的样本数N，可得到0到N－1的被赋予固定的频率间隔(等于频率分辨率(f0))下的振幅的值的功率谱。将这些功率谱中与最大值的功率谱对应的频率称为峰值频率。最大值的功率谱能通过公知的方法确定，例如，如果将1到N/2的功率谱逐个进行比较则能够确定。此时，如果在p点功率谱最大，则峰值频率fpk为fpk＝p×f0。
在专利文献1中，对如下方法进行了说明：同时满足频率分辨率小于等于12(Hz)、进行FFT的信号波的截取时间小于等于10(ms)(相当于水中的位置分辨率7.5m)，并求出峰值频率(参照该专利文献1的[0023]～[0026]、[0089]～[0090]。需要说明的是，在[0097]之后的数值例中，从输入信号波的截取时间变为5(ms)，但理由不明。)。当以频率分辨率12(Hz)优先时，根据[算式2]，截取时间为1/12(Hz)＝83.3(ms)。另一方面，在应用FFT的情况下，由于要求其样本数N是2的乘方，因此，在该专利文献1中，将信号截取时间设为相当于N＝1024的102.4(ms)(参照该专利文献1的[0055])。为了将进行FFT的输入信号波的截取时间设为102.4(ms)，由于以获得实际的数据的时间5(ms)是不够的，因此对不够的量102.4－5＝97.4(ms)期间的部分附加零值数据，将加入该零值数据所得的数据串(时间长度102.4(ms)的部分)进行FFT，求出峰值频率(参照该专利文献1的[0103])。
在该专利文献1记载的方法中，由于进行FFT的数据串中所含的有效数据串只不过是整体的约1/10(根据该专利文献1的[0103]，约1/20)，因此，当使施加于在输入信号波是模拟信号波时通常应用的A/D转换而得到的数字信号波的数字带通滤波器(在该专利文献1的图7中，数字BPF62)提高阶数并且为窄带时，通过数字带通滤波器后的信号迟钝，而且有效数据数量减少，因此必须减缓数字带通滤波器的滤波特性。就是说，在该专利文献1中，由于丝毫没有认识到难以避免时间长度短的输入信号波的情况下的由干扰噪声引起的负面影响这一潜在问题，因此，对那样的潜在问题的解决方案也没有任何启示和教导。
此外，在该专利文献1的情况下，在进行FFT的数字信号波(数字数据串)中，由于仅在应用FFT时所要求的时间窗长度约1/10(或者，约1/20)的长度的部分存在实际的数据，因此不得不产生大量的功率谱。在应用FFT的实际的信号波中包含干扰噪声的情况下，由于在峰值频率的周边进一步出现多余的频谱，因此存在难以确定峰值频率的可能性。
如此，当求取输入信号的峰值频率时，为了避免互反关系的制约而进行了各种尝试，但尚未找到常规的解决方法。
现有技术文献
专利文献
专利文献1：日本特开2012-247302号公报
发明内容
发明所要解决的问题
本发明的目的在于提供一种装置、方法以及程序，所述装置、方法以及程序能避免通过FFT进行频率分析时成为制约的、存在于频率分辨率(f0)与时间窗长度(T0)之间的f0＝1/T0这一互反关系，能以所希望的频率分辨率和所希望的时间窗长度来检测信号波的峰值频率。
用于解决问题的方案
(1)用于达成上述目的的峰值频率检测装置是检测规定的频带(fcl～fch)中功率谱最大的峰值频率的峰值频率检测装置，具备：n次方部，对数字数据串的各元素取n次方(n是大于等于2的整数)；FFT部，将与对根据采样频率fs、频率分辨率ftg、时间窗长度Ttg确定的N(N是2的乘方的整数)个取n次方后的采样频率fs的数字数据串进行快速傅立叶转换所得的功率谱的最大值对应的频率作为虚拟峰值频率导出；以及1/n倍部，将对所述虚拟峰值频率取1/n倍后的值作为数字数据串的峰值频率输出，所述峰值频率检测装置满足：
n≥1/(ftg×Ttg)
fs/(n×ftg)≤N≤fs×Ttg
fs＞2×n×fch。
需要说明的是，在本说明书中，fcl～fch是指大于等于fcl小于等于fch。
当假定作为频率分析的对象的信号波y是具有单一频率f(Hz)的[算式3]所示的正弦函数时，y的n次方(n是大于等于2的正整数)以[算式4]以及[算式5]表示。在[算式4]中n是奇数，在[算式5]中n是偶数。
[算式3]
y＝sin(2πft)
[算式4]
[算式5]
当对y取n次方(n＝2，3，4，……)时，无论n是奇数还是偶数，[算式4]以及[算式5]都表示：出现频率成分n×f，(n－2)×f，(n－4)×f，……。
因此，可以理解的是：如果能对y＝sin(2πft)取n次方，并从此取n次方得到的信号波中提取频率成分(n×f)的信号波，求出此信号波的峰值频率，那么如果对求得的频率取1/n倍，则其与原信号波的频率f对应。图1示意性地表示了该关系。
在此，将使用于频率分析的FFT的采样频率设为fs，将样本数设为N。当将y＝sin(2πft)进行该FFT并求取峰值频率时，则该峰值频率是f的计算值，此时的频率分辨率f0为：
[算式6]
f0＝fs/N。
在此，当将频率成分是(n×f)的正弦函数yn＝sin(2πnft)进行具有与求出y的峰值频率时相同值的fs和N的FFT并求取峰值频率时，求得的yn的频率fn’为(n×f)的计算值。在此，应该注意的是，该yn的频率fn’是在与关于y的频率分辨率相同的频率分辨率f0的基础上求得的。就是说，fn’的频率分辨率也是f0。因此，对fn’取1/n倍后的fn’/n是以频率分辨率f0/n求出y的频率f。
因此，当将直接进行FFT求出y的频率f时的频率分辨率设为f0，将时间窗长度设为T0，将根据yn求出y的频率f时的频率分辨率设为fn，将时间窗长度设为Tn时，则如下关系成立：
[算式7]
fn＝f0/n(n是大于等于2的整数)
[算式8]
Tn＝T0
[算式9]
fn×Tn＝(f0/n)×T0＝1/n＜1(n是大于等于2的整数)。
关于[算式9]，在图2中示出了n＝2，3，4时和不进行乘方时的f0×T0＝1的曲线图。
当求取原信号波的峰值频率时，频率分辨率ftg是用户所希望的峰值频率的频率分辨率。时间窗长度Ttg是用户所希望的FFT的时间窗长度。ftg与Ttg可以独立设定。
就是说，本发明是以满足如下条件的方式对原信号波的峰值频率进行计算的方法。
·峰值频率的频率分辨率fn≤ftg
·FFT的时间窗长度Tn≤Ttg
为此，根据[算式9]，需要满足：
[算式10]
fn×Tn＝1/n≤ftg×Ttg。
就是说，次方数n必须满足：
[算式11]
n≥1/(ftg×Ttg)。
需要说明的是，n不限于2的乘方2，4，8，16，32，……，只要选自大于等于2的整数即可。因此，可以根据用途来选择需要且足够的最优值。即，根据本发明，能避免通过FFT进行频率分析时成为制约的、存在于频率分辨率(f0)与时间窗长度(T0)之间的f0＝1/T0这一互反关系，能以所希望的频率分辨率和所希望的时间窗长度来检测信号波的峰值频率。具体而言，即使在ftg×Ttg＜1的范围，也能以所希望的频率分辨率和所希望的时间窗长度来检测信号波的峰值频率。
如[算式9]所示，将对信号波取n次方后的数字数据串进行FFT的情况下的峰值频率的频率分辨率fn、时间窗长度Tn有fn×Tn＝1/n的关系，如果对此进行图示，则如图4所示。该图4中所描绘的fn×Tn＝1/n的曲线图中的粗线部分为满足ftg、Ttg要求的fn、Tn。
当将fn＝ftg时的Tn设为Tmin时，则Tn的可取范围为：
[算式12]
Tmin≤Tn≤Ttg。
另一方面，由于
[算式13]
ftg×Tmin＝1/n
Tn＝N/fs，
因此，
[算式14]
1/(n×ftg)≤N/fs≤Ttg。
在进行FFT的情况下，由于N必须是2的乘方，因此N的可取范围为：
fs/(n×ftg)≤N≤fs×Ttg(N是2的乘方)。
假如在N不存在的情况下，增大n或增大fs、或者增大两者，直至N存在。不过，根据FFT的采样定理，必须满足：
[算式15]
fs＞2×n×fch。
在相应的N存在多个的情况下，无论采用哪个N都满足所希望的频率分辨率ftg和所希望的时间窗长度Ttg。不过，N越小，FFT的计算量越少，因此通常最好采用最接近下式的N。
[算式16]
fs/(n×ftg)
如此，即使在ftg×Ttg＜1的范围，也能确定能以所希望的频率分辨率ftg、时间窗长度Ttg来检测峰值频率的次方数n、FFT的采样频率fs、FFT的样本数N。需要说明的是，次方数n、fs、N的确定方法可以不按照上述过程。可以适当地应用数值反复进行试错直至满足条件，或者制作简单的程序来导出，使用什么样的方法都可以。
(2)用于达成上述目的的峰值频率检测装置具备：第一数字带通滤波器，提取所述规定的频带中所含的频率的数字数据串，所述第一数字带通滤波器的输出可以输入至所述n次方部。
当采用该构成时，由于能够为了从包含各种频率成分的信号中提取规定的频带的信号而简化或省略前一阶段的模拟滤波器，因此能缩小电路规模。
(3)用于达成上述目的的峰值频率检测装置具备：间隔剔除部，将采样频率fis的数字数据串间隔剔除为1/r(r是大于等于2的整数)，并将采样频率设为fs，所述间隔剔除部的输出可以输入至所述第一数字带通滤波器。
当采用该构成时，即使在输入至峰值频率检测装置的数字数据串的采样频率大于fs的情况下，也能以所希望的频率分辨率ftg、时间窗长度Ttg来检测峰值频率。
(4)用于达成上述目的的峰值频率检测装置具备：内插部，将接收数字数据串内插为g倍(g是大于等于2的整数)，并将采样频率设为fs，所述内插部的输出可以输入至所述第一数字带通滤波器。
当采用该构成时，即使在输入至峰值频率检测装置的数字数据串的采样频率小于fs的情况下，也能以所希望的频率分辨率ftg、时间窗长度Ttg来检测峰值频率。
(5)用于达成上述目的的峰值频率检测装置具备：第二数字带通滤波器，从取n次方后的N个所述数字数据串中提取第二频带中所含的数字数据串，由所述第二数字带通滤波器提取的数字数据串被输入至所述FFT部，所述第二频带可以是大致n×fcl～n×fch。
假设即使是仅具有单一频率成分的数字数据串，当取n次方时，也具有多个频率成分，因此在功率谱出现多个峰值，但所要检测的峰值出现在与频带fcl～fch对应的n×fcl～n×fch。因此，通过提取n×fcl～n×fch的成分，能检测规定的频带(fcl～fch)中功率谱最大的峰值频率。需要说明的是，在能从多个峰值频率中确定并选择所要检测的峰值频率的情况下，也可以不使用第二数字带通滤波器。
(6)用于达成上述目的的峰值频率检测装置具备：间隔剔除部，将由所述第二数字带通滤波器提取的数字数据串间隔剔除为1/r(r是大于等于2的整数)，并将采样频率设为fs，所述间隔剔除部的输出可以输入至所述FFT部。
当采用该构成时，即使在输入至峰值频率检测装置的数字数据串的采样频率大于fs的情况下，也能以所希望的频率分辨率ftg、时间窗长度Ttg来检测峰值频率。
(7)在用于达成上述目的的峰值频率检测装置中，具备：第一数字带通滤波器，提取所述特定的频带中所含的频率的数字数据串；第二数字带通滤波器，从所述n次方部的输出中提取第二频带中所含的数字数据串，所述第一数字带通滤波器的输出被输入至所述n次方部，所述第二数字带通滤波器的输出被输入至所述FFT部，所述第二频带可以是大致n×fcl～n×fch。
(8)用于达成上述目的的峰值频率检测装置，代替所述n次方部，具备多重乘方部，所述多重乘方部具备k(k是大于等于2的整数)级乘方块(j)，所述乘方块(j)具备：乘方部(j)(j＝1，2，……，k)，对所输入的数字数据串取mj次方；以及数字带通滤波器(j)，从所述乘方部(j)的输出中提取特定的频带fcl(j)～fch(j)的信号，可以是：
n＝m1×m2×……×mk
fcl(j)≈(m1×m2×……×mj)×fcl
fch(j)≈m1×m2×……×mj)×fch。
在对数字数据串进行乘方的情况下，次方数越小，越能减少由乘方产生的多余的频率成分。多余的频率成分越少，越容易通过数字带通滤波器来去除多余的频率成分。即，与能以一级的乘方部和第二数字带通滤波器应对的次方数n的上限相比，能以将乘方部和数字带通滤波器单元化并采用多级构成的方式应对的次方数n的上限更大。因此，通过采用该构成，能扩大可设定为次方数的n的范围。
(9)用于达成上述目的的峰值频率检测装置可以具备：操作部，受理用户的指示；以及参数设定部，对根据所述指示的n、fs、N的至少任一个进行设定。
通过采用该构成，用户无需以重复进行试错的方式设定满足下式的n、fs、N：
n≥1/(ftg×Ttg)
fs/(n×ftg)≤N≤fs×Ttg
fs＞2×n×fch。
需要说明的是，权利要求所述各部的功能通过功能由构成自身确定的硬件资源、功能由程序确定的硬件资源、或它们的组合来实现。此外，这些各部的功能并不限于通过各自物理上相互独立的硬件资源来实现。而且，本发明作为方法、作为计算机程序、作为计算机程序的记录介质都成立。当然，此计算机程序的记录介质既可以是磁记录介质，也可以是光磁记录介质，还可以是今后开发的任何记录介质。
附图说明
图1是表示本发明的实施方式的频谱的图。
图2是本发明的实施方式的曲线图。
图3是本发明的实施方式的框图。
图4是本发明的实施方式的曲线图。
图5是表示本发明的实施方式的频谱的图。
图6是本发明的实施方式的框图。
图7是本发明的实施方式的波形图。
图8是本发明的实施方式的框图。
图9是本发明的实施方式的框图。
图10是本发明的实施方式的框图。
图11是本发明的实施方式的波形图。
图12是本发明的实施方式的波形图。
图13是本发明的实施方式的框图。
图14是本发明的实施方式的框图。
图15是本发明的实施方式的框图。
图16是本发明的实施方式的画面构成图。
具体实施方式
以下，参照附图对本发明的实施方式进行说明。需要说明的是，将在各图中对应的构成要素赋予相同的附图标记，并省略重复的说明。
1、第一实施方式
在第一实施方式中，对峰值频率检测装置和使用了该装置的多普勒测量仪进行说明，所述峰值频率检测装置以满足所希望的频率分辨率ftg和所希望的时间窗长度Ttg的方式对以采样频率fs进行了采样的接收数字数据串的、由下限值fcl和上限值fch确定的频带fcl～fch中的峰值频率f进行检测。
如图3所示，作为本发明的第一实施方式的峰值频率检测装置1具备：第一数字带通滤波器(BPF)部11、n次方部12、第二数字带通滤波器(BPF)部13、FFT部14、以及1/n倍部15。
首先，对次方数n、使用于频率分析的FFT的采样频率fs、使用于频率分析的FFT的采样数N的确定方法的过程进行说明。
(步骤1：确定次方数n)
首先，确定满足n≥1/(ftg×Ttg)的条件的n(n是大于等于2的整数)。例如，将n设为满足n≥1/(ftg×Ttg)的最小的整数，如果在后述的步骤3中N不存在，则可以使n增加1并重新计算，也可以从最初就设得较大。
(步骤2：选择FFT的采样频率fs)
为了满足采样定理，以fs＞2×n×fch的方式来选择FFT的采样频率fs。
需要说明的是，在本实施方式中，输入至峰值频率检测装置1的数字数据串的采样频率为FFT的采样频率fs。在fs因电路上的制约等而不满足上式的情况下，使用其他实施方式。
(步骤3：确定FFT的样本数N)
接着，选择满足fs/(n×ftg)≤N≤fs×Ttg(N是2的乘方)的N。假如在N不存在的情况下，增大n或增大fs、或者增大两者，直至N存在。不过，根据FFT的采样定理，必须满足：
fs＞2×n×fch。
通过以上所说明的方法选择的次方数n、fs、N设定于n次方部12以及FFT部14。此外，要检测峰值频率的频带的下限值fcl和上限值fch被设定为第一数字BPF11的截止频率。此外，对于第二数字BPF13的截止频率，设定要检测峰值频率的频带的下限值fcl和上限值fch的n倍的值。
(根据数字数据串计算峰值频率)
当向设定了次方数n、fs、N的峰值频率检测装置1输入作为对象的数字数据串时，如下所述，以满足所希望的频率分辨率ftg和所希望的时间窗长度Ttg的方式检测出峰值频率。
当作为对象的数字数据串输入至峰值频率检测装置1时，第一数字BPF11从数字数据串的大范围的频率成分中排除所设定的频带以外的多余的直流成分、低频成分、高频成分，提取接近单一频率f的频率成分。
在作为要检测峰值频率的频带，fcl～fch过宽，其间存在多个大的功率谱的情况下，如图5所示，较好的是：重新将fcl～fch设定得较窄，分数次求出峰值频率。在图5的例子中，fcl～fch表示最初的频带，f1、f2表示具有大功率谱的两个频率，fcl’～fch’表示变窄后的频带。需要说明的是，fcl’～fch’的间隔无需是固定的，如图5所示，频率越高，越可以扩大fcl’～fch’的间隔。
在此，将从第一数字BPF11通过后的数字数据串设为A(1)：a0，a1，a2，…。在A(1)中，除了频率f，还包含功率谱小的一些多余的频率成分。
接着，将该A(1)输入n次方部12。在n次方部12中，将上述确定的n作为次方数对A(1)的各元素进行乘方。当将从n次方部12通过后的数字数据串设为B(n)：b0，b1，b2，……时，则bi＝(ai)n(i＝0，1，2，3，4，……)。在B(n)中，还包含：频率成分n×f，(n－2)×f，(n－4)×f，……这样的功率谱小的多余的低频成分和高频成分。
接着，将B(n)输入第二数字BPF13，提取n×fcl～n×fch的频带的频率成分。通过第二数字BPF13，能从B(n)的大范围的频率成分中排除多余的低频、高频成分，并提取接近单一频率的频率成分。虽然第二数字BPF13的频带设成第一数字BPF11的频带的n倍的n×fcl～n×fch为好，但根据用途，稍微改变也没有关系。将从第二数字BPF13通过后的数字数据串设为C(n)：c0，c1，c2，……。在C(n)中，除了频率n×f，还包含一些多余的频率成分。
在此，将输入至峰值频率检测装置1的数字数据串的采样频率fs设为fs＞2×n×fch的理由是C(n)的频带的上限是n×fch。就是说，接收数字数据串的采样频率fs以满足FFT的采样定理为条件。
接着，将C(n)输入FFT部14并计算出峰值频率。在FFT部14中，以所设定的采样频率fs和样本数N来对C(n)执行FFT，计算出峰值频率。将从FFT部14输出的峰值频率设为(n×f’)。
接着，使从FFT部14输出的峰值频率(n×f’)通过1/n倍部15，求出f’。求得的f’为输入至峰值频率检测装置1的数字数据串的峰值频率f的计算值。该f’是以满足ftg、Ttg的方式计算出的。
如以上所说明，峰值频率检测装置1能根据接收数字数据串，以所希望的频率分辨率ftg和所希望的时间窗长度Ttg来检测峰值频率。然后，由于峰值频率检测装置1不进行近似计算，因此不存在计算出的峰值频率的精度劣化。因此，能高精度地检测峰值频率。
图6是表示组装了峰值频率检测装置1的多普勒测量仪2的框图。在多普勒测量仪2中，具备：收发器21，其兼具向介质中发送信号波的功能和接收来自介质中的对象物的反射波的功能。收发器21介由收发切换电路22连接于发送电路23的输出端和接收放大器24的输入端。发送电路23生成发送频率ftx的信号。在接收放大器24的输出端设有将接收信号变频为中间频率信号的调制器25，本地振荡频率floc的信号从本地振荡电路26供给至调制器25。调制器25的输出、即中间频率信号介由模拟滤波器27输入至模拟/数字转换器(A/D)28，从而以FFT所要求的采样频率被转换为数字信号。接着，该数字信号输入至峰值频率检测装置1并通过FFT进行频率分析，从峰值频率检测装置1输出峰值频率。
在本实施例中，对在以下条件下检测海水中的声波的峰值频率的方法进行说明。
海水中的声波传播速度C：1500m/s
发送频率ftx：120kHz
本地发信频率floc：137kHz
检测最大速度(水平方向)V：15m/s
检测速度精度(水平方向)V0：0.15m/s
检测对象物的位置的精度(位置分辨率)D0：7.5m
此外，设为声波在相对于水平方向的倾斜方向(θ＝60°)上进行信息收发。
在C＞＞V的情况下，多普勒频率fdop为：
[算式17]
由于检测最大速度V是15m/s，因此：
[算式18]
多普勒信号是120±1.2kHz的范围，所观测的频带的宽度Δfp为：
Δfp＝2×1200Hz＝2400Hz。
在此，中间频率fmid设为floc－ftx＝137－120＝17kHz。
模拟滤波器27具有如下特性：
使fmid±(△fp/2)＝17000±1200Hz的信号通过，在下一阶段的A/D转换器28不产生混叠(aliasing)。
由于检测速度精度V0是0.15m/s，因此：
[算式19]
频率分辨率f0为12Hz。
由于检测对象物的位置的精度(位置分辨率)是7.5m，因此时间窗长度是作为往复7.5m的时间的、7.5×2/1500＝10ms。
即使直接将模拟/数字转换器(A/D)28的输出进行FFT，也不能满足该条件。其原因是，在频率分辨率是12Hz时，时间窗长度为1/12＝83.3ms(＞10ms)，位置分辨率为0.0833×1500/2＝62.5m(＞7.5m)。相反，当将时间窗长度设为10ms时，频率分辨率为1/0.01＝100Hz(＞12Hz)。
因此，通过具备n次方部12和1/n倍部15的峰值频率检测装置1来求出峰值频率。当将所希望的频率分辨率ftg设为12Hz，将所希望的时间窗长度Ttg设为10ms时，次方数n为：
n≥1/(ftg×Ttg)＝1/(12×0.01)＝8.3。在此，设为n＝12。
接着，确定A/D转换器28的采样频率fs。当第一数字BPF11的频带fcl～fch匹配17000±1200Hz时，fcl＝15800Hz，fch＝18200Hz。因此，A/D转换器28的采样频率fs为：
fs＞2×n×fch＝2×12×18200＝436800Hz，所以设为fs＝510kHz。
根据fs/(n×ftg)≤N≤fs×Ttg，FFT的样本数N是满足510000Hz/(12×12Hz)＝3541.7≤N≤510000Hz×0.01s＝5100的2的乘方的整数。即，N＝4096。
以这种方式确定采样频率fs和样本数N，将从模拟/数字转换器(A/D)28输出的数字数据串输入峰值频率检测装置1。
需要说明的是，作为第一数字BPF11，例如，可以优选使用阶数为8阶且将截止频率设定为fcl＝15.8kHz、fch＝18.2kHz的巴特沃斯型的IIR(Infinite impulse response：无限脉冲响应)滤波器。此外，作为第二数字BPF13，例如，可以优选使用阶数为8阶且将截止频率设定为189.6kHz(12×fcl)以及218.4kHz(12×fch)的巴特沃斯型的IIR滤波器。
在此，将对象物的相对移动速度(水平方向)V设为V＝10m/s，设为在相对于水平方向的倾斜方向(α＝60°)上对声波进行收发。
接收信号波的多普勒频率fdop为：
[算式20]
当将Sin{2π(17000+800)t}＝Sin(2π17800t)假定为向采样频率510kHz的A/D转换器28的输入信号，制成伪数字数据串后，通过上述设定实际求取峰值频率时，得到(12×f’)≈213662.1Hz。所求的f’为：f’＝(12×f’)/12＝213662.1/12≈17805.2Hz。f’是以满足ftg的方式计算出的。实际上，对该数字数据串的测量误差ε是ε＝f’－fdop＝17805.2－17800＝5.2Hz，包括在±(ftg)/2＝±6Hz以下。需要说明的是，由于频率分辨率f0时的峰值频率为f0的间隔中最近的点，因此频率分辨率f0时的峰值频率的误差为±(f0/2)以下。
在本实施例中，将从第二数字BPF13输出的数字数据串的实例在图7中示出。在图7中，T1是进行FFT的区间。T2是直至数字数据串的振幅稳定的区间。由于T1的元素数是4096个，因此T1的区间长度为4096/510000≈8.0ms，满足所希望的时间窗长度Ttg＝10ms。
需要说明的是，更优选将区间T2中振幅稳定后的数字数据串进行FFT。由于T2的元素数是约600个，因此T2的区间长度为600/510000≈1.2ms。即使是在这种情况下，由于将次方数n设定为较大的12，因此计算峰值频率所需的数字数据串的时间长度是1.2+8.0＝9.2ms，小于等于Ttg＝10ms。需要说明的是，在本实施例中，即使将T2之后的4096个数字数据串进行FFT，计算结果也是f’≈17805.2Hz，得到相同值。
如此，通过增大次方数n，能轻松地满足Ttg，增大第一数字BPF11、第二数字BPF13的阶数使其陡峭，能抵抗外来噪声。需要说明的是，可容易地理解的是：即使Ttg是5ms，通过增大次方数n，也能以所希望的频率分辨率ftg来进行频率分析。
在此，对设置作为带通滤波器的第一数字BPF11的理由进行说明。为了应用[算式4]以及[算式5]，在取n次方之前，需要去除频带fcl～fch以外的多余的直流成分、低频成分、高频成分。如果能通过模拟滤波器27来实现，则可以不需要第一数字BPF11。但是，这样的模拟滤波器为高阶且高精度，电路规模大，成本高。因此，模拟滤波器27采用在A/D转换器28不产生混叠的程度的设计，频带fcl～fch的成分的提取采用能容易地实现高阶且高精度的设计并且成本低的数字带通滤波器才是上策。
2、第二实施方式
图8是表示作为本发明的第二实施方式的峰值频率检测装置3的构成的框图。峰值频率检测装置3为在峰值频率检测装置1的第一数字BPF11之前追加有间隔剔除部16的构成。这是为了通过对A/D转换后的数字数据串进行间隔剔除使其减少来降低采样频率。在输入至峰值频率检测装置3的数字数据串的采样频率fs高、次方数n和第一数字BPF11的截止频率fch满足fs＞4×n×fch的情况下，优选应用第二实施方式。
在间隔剔除部16中，将数字数据串间隔剔除为1/r(r是大于等于2的整数)，间隔剔除后的数字数据串为：
fs＞2×n×fch。fs是间隔剔除后的采样频率。
对于间隔剔除为1/r，例如，以如下方式进行：将间隔剔除前的数字数据串设为P(1)：p0，p1，p2，……，将间隔剔除后的数字数据串设为Q(1)：q0，q1，q2，……时，进行qi＝p(r×i)(i＝0，1，2，3，4，……)。在r＝2时，为q0＝p0，q1＝p2，q2＝p4，q3＝p6。
间隔剔除的方法可以是其他方法。例如，在r＝2时，也可以采用q0＝(p0+p1)/2，q1＝(p2+p3)/2，q2＝(p4+p5)/2，……等。
需要说明的是，在设为r＝2时，不存在满足fs/(n×ftg)≤N≤fs×Ttg(N是2的乘方)的N，也无法增大n的情况下，无法应用本实施例。该情况应用第一实施方式。
以下，对将本实施例的峰值频率检测装置3应用于图6所示的多普勒测量仪2的情况进行说明。需要说明的是，下述数值数据以外所需的数值数据设为与第一实施方式的数值例相同。此外，为了简化之后的说明，与第一实施方式相同，当将fs＝510kHz的A/D转换器28的采样频率设为10.2MHz时，输入至间隔剔除部16的接收数字数据串的采样频率fis为10.2MHz。
在此，fis＝10.2MHz＞4×n×fch＝4×12×18.2kHz＝873.6kHz。另一方面，2×n×fch＝2×12×18.2kHz＝436.8kHz。因此，设为将接收数字数据间隔剔除为1/r＝1/20。当将间隔剔除后的数字数据串的采样频率设为fs时，
fs＝fis/r＝10.2MHz/20＝510kHz＞436.8kHz。
当将间隔剔除前的数字数据串设为P(1)：p0，p1，p2，……，将间隔剔除后的数字数据串设为Q(1)：q0，q1，q2，……时，在间隔剔除部16中，进行q0＝p0，q1＝p20，q2＝p40，q3＝p60，……。
可以理解的是：该Q(1)与以采样频率fs＝510kHz对输入信号进行A/D转换后的信号串同等。因此，FFT的采样频率fs为fs＝510kHz。
与第一实施方式相同，在将y＝sin(2π17800t)假定为向采样频率10.2MHz的A/D转换器28的输入信号，制成伪数据串后，以上述实施例的设定实际进行频率分析的结果与第一实施方式的数值例相同，满足ftg、Ttg。
3、第三实施方式
图9是表示本发明的第三实施方式的峰值频率检测装置4的构成的框图。峰值频率检测装置4为在第一实施方式的第二数字BPF13之后追加有间隔剔除部17的构成。对于第一数字BPF部11、n次方部12、第二数字BPF13，只不过是所输入的数字数据的采样频率在第一实施方式中为fs而在此为fis，动作是相同的。在接收数字数据串的采样频率fs高、相对于次方数n和第一数字BPF11的截止频率fch，fs＞4×n×fch的情况下，为了削减计算量，一般应用第二实施方式为好。但是，即使是在本实施方式中，也能进行满足所希望的时间窗长度ftg、所希望的时间窗长度Ttg的处理。
在间隔剔除部17中，在将接收数字数据串间隔剔除为1/r(r是大于等于2的整数)，将间隔剔除后的数字数据串的采样频率设为fs时，
fs＞2×n×fch。
间隔剔除后的数字数据串输入至FFT部14，此后的处理与第一实施方式相同。
需要说明的是，在不满足r≥2且fs/(n×ftg)≤N≤fs×Ttg(N是2的乘方)的情况下，无法应用本实施方式。该情况下，应用第一实施方式。
以下，对将本实施例的峰值频率检测装置4应用于图6所示的多普勒测量仪2的情况进行说明。需要说明的是，下述数值数据以外所需的数值数据设为与第一实施方式的数值例相同。在A/D转换器28的采样频率是10.2MHz的情况下，输入至间隔剔除部17的接收数字数据串的采样频率fis为10.2MHz。
在此，fis＝10.2MHz＞4×n×fch＝4×12×18.2kHz＝873.6kHz。另一方面，2×n×fch＝2×12×18.2kHz＝436.8kHz。因此，设为将从第二数字BPF13通过后的数字数据间隔剔除为1/r＝1/20。当将间隔剔除后的数字数据串的采样频率设为fs时，fs＝fis/r＝10.2MHz/20＝510kHz＞436.8kHz。FFT的采样频率fs为fs＝510kHz。
与第一实施方式的数值例相同，在将y＝sin(2π17800t)假定为向采样频率10.2MHz的A/D转换器28的输入信号，制成伪数字数据串后，以上述实施例的设定实际进行频率分析的结果与第一实施方式的数值例相同，满足ftg、Ttg。
4、第四实施方式
(以下，在申请时将段落号上调)图10是表示作为本发明的第四实施方式的峰值频率检测装置5的构成的框图。峰值频率检测装置5为在第一实施方式的数字BPF11之前追加有内插部18的构成。通过追加内插部18，能对A/D转换后的数字数据串进行内插使其增加，由此提高采样频率。在输入至峰值频率检测装置5的数字数据串的采样频率fs低、相对于次方数n和第一数字BPF11的较高的截止频率fch，2×fch＜fs＜2×n×fch的情况下，应用本实施方式。
在内插部18中，将接收数字数据串内插为g倍(g是大于等于2的整数)，从而内插后的数字数据串的采样频率fs为fs＞2×n×fch。例如，在将内插前的数字数据串设为U(1)：u0，u1，u2，……，将内插后的数字数据串设为V(1)：v0，v1，v2，……时，内插部18以如下方式进行内插：
vi＝u0(i＝0，1，2，3，……，(g－1))
vi＝u1(i＝g，g+1，g+2，g+3，……，(2g－1))
vi＝u2(i＝2g，2g+1，2g+2，2g+3，……，(3g－1))
……
在g＝2时，v0＝u0，v1＝u0，v2＝u1，v3＝u1，v4＝u2，v5＝u2，……。内插的方法可以是其他方法。例如，在g＝2时，也可以以如下方式进行内插：
v0＝u0
v1＝(u0+u1)/2
v2＝u1
v3＝(u1+u2)/2
v4＝u2
v5＝(u2+u3)/2
……
以下，对将本实施例的峰值频率检测装置5应用于图6所示的多普勒测量仪2的情况进行说明。需要说明的是，下述数值数据以外所需的数值数据设为与第一实施方式的数值例相同。为了简化说明，与第一实施方式相同，将fs＝510kHz的A/D转换器28的采样频率设为42.5kHz。此时，输入至内插部18的接收数字数据串的采样频率fis为42.5kHz。该情况下，由于：
2×fch＝2×18.2kHz＝36.4kHz
2×n×fch＝2×12×18.2kHz＝436.8kHz
2×fch＜fis＜2×n×fch
因此应用本实施方式。
为了将接收数字数据串的采样频率扩大为大于等于436.8kHz，设成内插为g＝12倍。当将内插后的数字数据串的采样频率设为fs时，
fs＝fis×g＝42.5kHz×12＝510kHz＞436.8kHz。
在将内插前的数字数据串设为U(1)：u0，u1，u2，……，将内插后的数字数据串设为V(1)：v0，v1，v2，……时，内插部18以如下方式进行内插：
v0～v11＝u0
v12～v23＝u1
v24～v35＝u2
v36～v47＝u3
……
当对内插后的数字数据串V(1)的例子进行图示时，如图11所示，呈阶梯状。图11是将y＝sin(2π17000t)假定为向采样频率42.5kHz的A/D转换器28的输入信号，在制成伪数字数据串后进行内插的例子。该V(1)的采样频率fs为510kHz。因此，此后的处理与图10的构成的情况相同，将V(1)视为向第一数字BPF11输入的采样频率fs的接收数字数据串，只要以相同的方法进行频率分析即可。
与第一实施方式的数值例相同，将y＝sin(2π17800t)假定为向采样频率42.5kHz的A/D转换器28的输入信号，在制成伪数字数据串后，以上述实施例的设定实际进行频率分析的结果与第一实施方式的数值例相同，满足ftg、Ttg。
在此，参照图12，对将n次方部12之前的数字滤波器设为作为带通滤波器的数字带通滤波器的理由进行说明。
图12A是向A/D转换器28的输入信号，设为：
y＝sin(2π17000t)+2。
就是说，设为：残留有未被模拟滤波器滤净的直流成分2的、振幅1、频率17kHz的sin波形。以采样频率为17kHz的2.5倍的42.5kHz对该输入信号进行A/D转换，其相当于输入至峰值频率检测装置1的数字数据串。
通过内插部18将该数字数据串扩大为12倍。图12B是从内插部18输出的数字数据串的例子。该数字数据串的采样频率变为相当于12×42.5kHz＝510kHz。
图12C是在使用了阶数为8阶且将截止频率设定为fcl＝15.8kHz、fch＝18.2kHz的巴特沃斯型的IIR滤波器的情况下，从第一数字BPF部11输出的数字数据串。第一数字BPF部11的输出被消除了从内插部18输出的多余的直流成分、多余的低频成分、多余的高频成分，收敛后的波形成为接近sin波形的波形。
如此，如果预先通过第一数字BPF部11来限制频带，则在接下来的n次方部12取n次方时，能抑制[算式4]、[算式5]以外的频率成分的产生。即，在n次方部12的前一阶段具备数字带通滤波器，由此即使数字数据串含有多余的频率成分也能应对。因此，能使前一阶段的模拟滤波器27的性能要求降低至在下一阶段的A/D转换器28不产生混叠的程度，也能使A/D转换器28的采样频率降低，因此能削减电路的规模和成本。
需要说明的是，截止频率fcl、fch的数字带通滤波器也可以采用将截止频率fcl的数字高通滤波器与截止频率fch的数字低通滤波器进行组合的构成。
5、第五实施方式
图13是表示作为本发明的第五实施方式的峰值频率检测装置6的构成的框图。峰值频率检测装置6是将第一实施方式的n次方部12和第二数字BPF13替换为多重乘方部19的装置。
如图14所示，多重乘方部19按照k级(k是大于等于2的整数)序号顺序级联有乘方块(j)，该乘方块(j)包括：乘方部(j)(j＝1，2，3，……，k)，对数字数据串取mj次方(mj是大于等于2的整数)；以及数字带通滤波器(j)(BPF(j))，从乘方部(j)的输出中提取特定的频带fcl(j)～fch(j)的信号。
乘方块(j)的乘方部(j)的次方数mj以如下方式进行选择：关于替换前的n次方部的次方数n，
n＝m1×m2×……×mk成立。
此外，设定为：
fcl(j)≈(m1×m2×……×mj)×fcl
fch(j)≈(m1×m2×……×mj)×fch
在次方数mj是偶数的情况下，根据[算式5]，当对y取mj次方时，产生直流成分。因此，需要数字BPF(j)能消除直流成分。
根据[算式4]、[算式5]可知：在对sin(2πft)取m次方的情况下，如果减小次方数m，则能减少所产生的小于等于(m－2)×f的频率成分的数量。因此，更容易通过其后的数字BPF来抑制所需频带以外的频率成分。因此，当应用本实施方式时，能应对至更高的次方数n。作为目标，在n次方部的次方数n大于16的情况下，考虑采用本实施方式为好。
在本实施方式中，与第一实施方式相同，也满足所希望的频率分辨率ftg和所希望的时间窗长度Ttg。此外，在本实施方式中，未进行近似计算，因此也不存在所计算的峰值频率的精度劣化。
在本实施方式中，由于将乘方部和第二数字BPF部设为多级构成，因此与第一实施方式相比计算量增加。但是，能够满足所希望的频率分辨率ftg和所希望的时间窗长度Ttg并且无精度劣化地对峰值频率进行计算直至更高的次方数n这一优点弥补这些负面影响绰绰有余。
以下，对将本实施例的峰值频率检测装置6应用于图6所示的多普勒测量仪2的情况进行说明。需要说明的是，下述数值数据以外所需的数值数据设为与第一实施方式的数值例相同。为了简化之后的说明，与第一实施方式相同，设为：n＝12、fs＝510kHz、N＝4096、n＝m1×m2。此外，设为：m1＝4，m2＝3。
在将y＝sin(2π17800t)假定为向采样频率510kHz的A/D转换器28的输入信号并制成伪数字数据串后，向第一数字BPF11输入，将第一数字BPF11的输出设为A(1)：a0，a1，a2，……。
在乘方部(1)，对A(1)的各元素取4次方，将从乘方部(1)通过后的数字数据串设为B(4)：b0，b1，b2，……。即，bi＝(ai)4(i＝0，1，2，3，4，……)。
使该B(4)从数字BPF(1)通过，并将从数字BPF(1)数字通过后的数字数据串设为C(4)：c0，c1，c2，……。数字BPF(1)设为阶数为8阶且截止频率63.2kHz(4×fcl)、72.8kHz(4×fch)的巴特沃斯型的IIR滤波器。
接着，在乘方部(2)，对C(4)的各元素取3次方，将从乘方部(2)通过后的数字数据串设为D(12)：d0，d1，d2，……。即，di＝(ci)3(i＝0，1，2，3，4，……)。
使该D(12)从数字BPF(2)通过，并将从数字BPF(2)通过后的数字数据串设为E(12)：e0，e1，e2，……。数字BPF(2)设为阶数为8阶且截止频率189.6kHz(4×3×fcl)、218.4kHz(4×3×fch)的巴特沃斯型的IIR滤波器。
关于该E(12)，当以样本频率fs＝510kHz、样本数N＝4096的FFT来求取频率时，得出(12×f’)≈213662.1Hz。所求的f’为：f’＝(12×f’)/12≈17805.2Hz。f’是以满足ftg的方式计算出的。实际上，在此的相对于理论值f＝17800Hz的误差ε是：ε＝f’－f＝17805.2－17800＝5.2Hz，包括在±(ftg)/2＝±12/2＝±6Hz。
此外，与第一实施方式相同，时间窗长度也是N/fs＝4096/510kHz＝8.0ms＜Ttg＝10ms，满足所希望的时间窗长度Ttg。
6、第六实施方式
可以在此前所说明的峰值频率检测装置1、3～6追加参数设定部。图15是表示在第一实施方式的峰值频率检测装置1追加了参数设定部20的构成的框图。参数设定部20是具备处理器、存储器、输入输出机构的计算机，根据使用了未图示的键盘、鼠标、触摸屏显示器等未图示的操作部的用户的输入，在峰值频率检测装置1、3～6设定参数的值。参数如此前所说明，如下：
接收数字数据串的采样频率fis
数字数据串的采样频率fs
所希望的频率分辨率ftg
所希望的时间窗长度Ttg
n次方部的次方数n(n是大于等于2的整数)
第一数字BPF的频带大致fcl～fch(fcl＜fch)
第二数字BPF的频带大致n×fcl～n×fch
FFT的采样频率fs
FFT的采样数N
参数设定部20既可以将所述全部参数的数值与用户可选择的输入值对应地预先存储在存储器中，并根据用户的输入值来进行设定，也可以将一部分参数的数值与用户可选择的输入值对应地预先存储在存储器中，并根据用户的输入来进行设定，剩余的参数根据所输入的数值由处理器通过计算求出。
在本实施方式中，对用作调音辅助装置的峰值频率检测装置进行说明，所述调音辅助装置使用户仅输入与频率对应的音阶序号，根据所输入的音阶序号从存储器取得全部参数的数值并进行设定。即，参数设定部20从存储器取得与音阶序号1～88中任一个的输入对应的参数的数值并进行设定。图16是输入音阶序号状态下的触摸屏显示器的画面构成图的一个例子。全部参数的数值事先根据音阶序号P全部确定并存储在非易失性存储器中。音阶序号与平均律钢琴的琴键对应，音阶序号1与最低音的频率对应，音阶序号88与最高音的频率对应。即，对于与音阶序号P对应的频率fp，由下述[算式21]确定的值被预先存储在非易失性存储器中。
[算式21]
然后，例如，相对于P＝49(fp＝442Hz)的各参数的数值以如下方式确定。
[算式22]
fs＝24kHz
Ttg＝＝(1/0.510913481)/10＝0.195727855s
n＝16
n×fc1＝16×fc1＝6870.673888Hz
n×fch＝16×fch＝7279·225418Hz
N＝4096
需要说明的是，ftg是相当于fp的2音分(cent)的量的频率。众所周知，音分值是将两个音之间的频率比以对数形式表示的值，100音分相当于平均律12音阶的半音。Ttg是时间窗长度(1/ftg)的1/10。fcl是fp的50音分下的频率。fch是fp的50音分上的频率。N必须是fs/(n×ftg)＝24000/(16×0.5109)≈2936～fs×Ttg＝24000×0.1957≈4697之间的2的乘方的整数。
如上所述，如果事先确定，参数设定部20就能根据音阶序号P的输入来设定全部参数的数值。通过使用如此设定的参数，能计算接收数字数据串的峰值频率。需要说明的是，向峰值频率检测装置1输入的数字数据串既可以通过未图示的话筒、A/D转换器实时地依次输入至峰值频率检测装置1，也可以存储于存储器。
在设为P＝49(fp＝442Hz)，并制成以采样频率24kHz对y＝sin(2π442t)进行采样的伪数字数据串后，实际进行频率分析的结果是得到16×f’≈7072.2656Hz。所求的f’为f’＝(16×f’)/16≈442.0166Hz。求得的峰值频率的频率分辨率满足所希望的频率分辨率ftg。实际上，在此的相对于理论值f＝442Hz的误差ε是ε＝f’－f＝442.0166－442＝0.0166Hz，包括在±(ftg/2)＝±0.5109/2≈±0.255Hz。此外，时间窗长度是N/fs＝4096/24kHz＝0.171s，满足所希望的时间窗长度Ttg＝0.195s。需要说明的是，求得的峰值频率的输出形式可以是单位Hz，也可以是相对于fp＝442Hz的音分值(+0.065音分)。
7、效果
根据以上所说明的本发明的实施方式，能避免通过FFT进行频率分析时成为制约的、存在于频率分辨率(f0)与时间窗长度(T0)之间的f0＝1/T0这一互反关系的问题，能以所希望的频率分辨率ftg和所希望的时间窗长度Ttg来检测信号波的峰值频率。如果接收数字数据串的采样频率大于等于2×fch，则能计算峰值频率。然后，由于无需进行近似计算、曲线拟合、取平均等处理，因此不存在峰值频率的计算精度的劣化。
此外，即使在输入信号包含多余的频带的成分(直流成分、低频成分，高频成分)、或者S/N比差、或者几乎没有振幅的情况下，由于通过第一数字BPF、第二数字BPF来提取计算峰值频率所需的频带的频率成分，因此能毫无问题地计算出峰值频率。因此，能降低峰值频率检测装置的前一阶段的硬件的要求规格，能谋求小型化、降低成本。具体而言，例如，相对于第一数字BPF的截止频率fch，如果输入信号的采样频率大于等于2×fch，则能计算峰值频率，因此能降低前一阶段的模拟滤波器、A/D转换器等硬件的要求规格。此外，如果设定于乘方部的次方数n设定得较大，则能提高数字带通滤波器的阶数使其陡峭，抵抗外来噪声。
8、其他实施方式
需要说明的是，本发明的技术范围并不限定于上述实施例，在不脱离本发明主旨的范围内可以进行各种变更。
例如，虽然举例示出了IIR、巴特沃斯型、8阶的数字BPF，但也可适当使用其他形式的数字BPF。例如，作为能优选使用的例子，可列举出FIR(Finite impulse response：有限脉冲响应)、切比雪夫型等。阶数也不限于8阶。此外，截止频率也是根据情况扩大或缩小提取范围即可。此外，如果输入信号是以0为中心振荡的、几乎没有谐波成分的信号(就是说，接近单一频率的sin曲线的信号)，则可以将第一数字BPF设为低通滤波器。
此外，以反射回波中的多普勒频率的检测、向调音辅助装置的应用为例对本发明进行了说明，但本发明的应用范围不限于此。一般情况下，本发明能广泛应用于通过FFT进行信号波的峰值频率的检测的用途。
此外，上述实施方式的各功能部可以通过一个或多个LSI(Large ScaleIntegration：大规模集成电路)来实现，此外，多个功能部可以通过一个LSI来实现。此外，作为集成化的方法并不限定于LSI，也可以通过专用电路或通用处理器来实现。而且，可以在LSI制造后，利用可编程的FPGA(Field Programmable Gate Array：现场可编程门阵列)、可重新构成LSI内部的电路元件的连接或设定的可重构处理器(ReconfigurableProcessor)。此外，如果通过半导体技术的进步或派生的其他技术，出现替代LSI的集成电路化的技术，当然也可以使用此技术来进行功能块的集成化。
附图标记说明：
1 峰值频率检测装置
2 多普勒测量仪
3 峰值频率检测装置
4 峰值频率检测装置
5 峰值频率检测装置
6 峰值频率检测装置
11 第一数字BPF部
12 n次方部
13 第二数字BPF部
14 FFT部
15 1/n倍部
16 间隔剔除部
17 间隔剔除部
18 内插部
19 多重乘方部
20 参数设定部
21 收发器
22 收发切换电路
23 发送电路
24 接收放大器
25 调制器
26 本地振荡电路
27 模拟滤波器
28 A/D转换器

Claims (11)

1.一种峰值频率检测装置，检测规定的频带(fcl～fch)中功率谱最大的峰值频率，具备：
n次方部，对数字数据串的各元素取n次方(n是大于等于2的整数)；
FFT部，将与对根据采样频率fs、频率分辨率ftg、时间窗长度Ttg确定的N(N是2的乘方的整数)个取n次方后的采样频率fs的数字数据串进行快速傅立叶转换所得的功率谱的最大值对应的频率作为虚拟峰值频率导出；以及
1/n倍部，将对所述虚拟峰值频率取1/n倍后的值作为数字数据串的峰值频率输出，
所述峰值频率检测装置满足：
n≥1/(ftg×Ttg)
fs/(n×ftg)≤N≤fs×Ttg
fs＞2×n×fch。
2.根据权利要求1所述的峰值频率检测装置，其中，
具备：第一数字带通滤波器，提取所述规定的频带中所含的频率的数字数据串，
所述第一数字带通滤波器的输出被输入至所述n次方部。
3.根据权利要求2所述的峰值频率检测装置，其特征在于，
具备：间隔剔除部，将采样频率fis的数字数据串间隔剔除为1/r(r是大于等于2的整数)，并将采样频率设为fs，
所述间隔剔除部的输出被输入至所述第一数字带通滤波器。
4.根据权利要求2所述的峰值频率检测装置，其特征在于，
具备：内插部，将所述数字数据串内插为g倍(g是大于等于2的整数)，并将采样频率设为fs，
所述内插部的输出被输入至所述第一数字带通滤波器。
5.根据权利要求1所述的峰值频率检测装置，
具备：第二数字带通滤波器，从取n次方后的N个所述数字数据串中提取第二频带中所含的数字数据串，
由所述第二数字带通滤波器提取的数字数据串被输入至所述FFT部，
所述第二频带是大致n×fcl～n×fch。
6.根据权利要求5所述的峰值频率检测装置，
具备：间隔剔除部，将由所述第二数字带通滤波器提取的数字数据串间隔剔除为1/r(r是大于等于2的整数)，并将采样频率设为fs，
所述间隔剔除部的输出被输入至所述FFT部。
7.根据权利要求1所述的峰值频率检测装置，具备：
第一数字带通滤波器，提取所述特定的频带中所含的频率的数字数据串；以及
第二数字带通滤波器，从所述n次方部的输出中提取第二频带中所含的数字数据串，
所述第一数字带通滤波器的输出被输入至所述n次方部，
所述第二数字带通滤波器的输出被输入至所述FFT部，
所述第二频带是大致n×fcl～n×fch。
8.根据权利要求1～4中任一项所述的峰值频率检测装置，
代替所述n次方部，具备多重乘方部，所述多重乘方部具备k(k是大于等于2的整数)级乘方块(j)，所述乘方块(j)具备：乘方部(j)(j＝1，2，……，k)，对所输入的数字数据串取mj次方(mj是大于等于2的整数)；以及数字带通滤波器(j)，从所述乘方部(j)的输出中提取特定的频带fcl(j)～fch(j)的信号，
n＝m1×m2×……×mk
fcl(j)≈(m1×m2×……×mj)×fcl
fch(j)≈m1×m2×……×mj)×fch。
9.根据权利要求1～8中任一项所述的峰值频率检测装置，具备：
操作部，受理用户的指示；以及
参数设定部，对根据所述指示的n、fs、N中的至少任一个进行设定。
10.一种峰值频率检测方法，检测规定的频带(fcl～fch)中功率谱最大的峰值频率，包含：
对数字数据串的各元素取n次方(n是大于等于2的整数)；
将与对根据采样频率fs、频率分辨率ftg、时间窗长度Ttg确定的N(N是2的乘方的整数)个取n次方后的采样频率fs的数字数据串进行快速傅立叶转换所得的功率谱的最大值对应的频率作为虚拟峰值频率导出；以及
将对所述虚拟峰值频率取1/n倍后的值作为数字数据串的峰值频率输出，
所述峰值频率检测方法满足：
n≥1/(ftg×Ttg)
fs/(n×ftg)≤N≤fs×Ttg
fs＞2×n×fch。
11.一种峰值频率检测程序，检测规定的频带(fcl～fch)中功率谱最大的峰值频率，使计算机发挥功能以作为：
n次方部，对数字数据串的各元素取n次方(n是大于等于2的整数)；
FFT部，将与对根据采样频率fs、频率分辨率ftg、时间窗长度Ttg确定的N(N是2的乘方的整数)个取n次方后的采样频率fs的数字数据串进行快速傅立叶转换所得的功率谱的最大值对应的频率作为虚拟峰值频率导出；以及
1/n倍部，将对所述虚拟峰值频率取1/n倍后的值作为数字数据串的峰值频率输出，所述峰值频率检测程序满足：
n≥1/(ftg×Ttg)
fs/(n×ftg)≤N≤fs×Ttg
fs＞2×n×fch。
CN201480079481.8A 2014-07-10 2014-07-10 峰值频率检测装置以及方法 Expired - Fee Related CN106415310B (zh)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2014/068432 WO2016006079A1 (ja) 2014-07-10 2014-07-10 ピーク周波数検出装置、方法およびプログラム

Publications (2)

Publication Number Publication Date
CN106415310A true CN106415310A (zh) 2017-02-15
CN106415310B CN106415310B (zh) 2018-02-02

Family

ID=53888036

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201480079481.8A Expired - Fee Related CN106415310B (zh) 2014-07-10 2014-07-10 峰值频率检测装置以及方法

Country Status (6)

Country Link
US (1) US9857399B2 (zh)
EP (1) EP3168638A4 (zh)
JP (1) JP5766896B1 (zh)
KR (1) KR101785714B1 (zh)
CN (1) CN106415310B (zh)
WO (1) WO2016006079A1 (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107589424A (zh) * 2017-08-31 2018-01-16 努比亚技术有限公司 超声波采样方法、装置及计算机可读存储介质
CN114371342A (zh) * 2022-03-21 2022-04-19 国仪量子(合肥)技术有限公司 Fpga及基于其的实时信号测频方法以及锁相放大器
CN118226123A (zh) * 2024-02-27 2024-06-21 青岛汉泰电子有限公司 一种用于示波器中快速查询fft峰值的方法

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107402326B (zh) * 2017-07-20 2019-08-23 南京理工大学 一种改进s变换的有限窗长时频分析方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1143863A (zh) * 1994-11-22 1997-02-26 三星电子株式会社 通过窗口滤波器的直接扩频通信系统的接收机
JPH09298572A (ja) * 1996-04-30 1997-11-18 Advantest Corp Ｐｓｋ信号の搬送波周波数推定器
CN1975756A (zh) * 2005-11-30 2007-06-06 U－Comm.技术有限公司 包括峰值检测装置的rfid系统
JP2009192516A (ja) * 2008-01-16 2009-08-27 Mitsubishi Electric Corp 伝搬遅延時間測定装置及びレーダ装置
JP2012247303A (ja) * 2011-05-27 2012-12-13 Sonic Corp 周波数検出方法及び装置

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH01110351A (ja) * 1987-10-23 1989-04-27 Aloka Co Ltd 超音波ドプラ診断装置
JP2808954B2 (ja) * 1991-11-13 1998-10-08 国際電信電話株式会社 無変調信号検出及び周波数引き込み装置
US5431169A (en) * 1992-08-03 1995-07-11 Olympus Optical Co., Ltd. Ultrasonic diagnosing apparatus
JP3271504B2 (ja) * 1996-02-02 2002-04-02 三菱電機株式会社 周波数推定回路およびそれを用いたａｆｃ回路
US5799038A (en) 1996-04-30 1998-08-25 Advantest Corporation Method for measuring modulation parameters of digital quadrature-modulated signal
US7239674B2 (en) * 2003-06-04 2007-07-03 Honeywell Federal Manufacturing & Technologies, Llc Method of differential-phase/absolute-amplitude QAM
US20080052335A1 (en) * 2006-08-01 2008-02-28 Gee Edward C Systems and methods for time domain to frequency domain conversion using frequency shifting
JP4352082B2 (ja) * 2007-06-18 2009-10-28 株式会社東芝 周波数同期回路、方法、プログラム及びこれらを用いた受信装置
US8390508B1 (en) 2010-04-05 2013-03-05 Raytheon Company Generating radar cross-section signatures
JP5094937B2 (ja) 2010-09-21 2012-12-12 三菱電機株式会社 周波数解析装置
JP5670836B2 (ja) 2011-05-27 2015-02-18 株式会社ソニック フーリエ変換でのサンプル数を削減した、短時間信号のピークパワースペクトルを検出する方法及び装置
JP2012247304A (ja) * 2011-05-27 2012-12-13 Sonic Corp 短時間信号のピークパワースペクトルを検出する方法及び装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1143863A (zh) * 1994-11-22 1997-02-26 三星电子株式会社 通过窗口滤波器的直接扩频通信系统的接收机
JPH09298572A (ja) * 1996-04-30 1997-11-18 Advantest Corp Ｐｓｋ信号の搬送波周波数推定器
CN1975756A (zh) * 2005-11-30 2007-06-06 U－Comm.技术有限公司 包括峰值检测装置的rfid系统
JP2009192516A (ja) * 2008-01-16 2009-08-27 Mitsubishi Electric Corp 伝搬遅延時間測定装置及びレーダ装置
JP2012247303A (ja) * 2011-05-27 2012-12-13 Sonic Corp 周波数検出方法及び装置

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107589424A (zh) * 2017-08-31 2018-01-16 努比亚技术有限公司 超声波采样方法、装置及计算机可读存储介质
CN107589424B (zh) * 2017-08-31 2021-12-17 努比亚技术有限公司 超声波采样方法、装置及计算机可读存储介质
CN114371342A (zh) * 2022-03-21 2022-04-19 国仪量子(合肥)技术有限公司 Fpga及基于其的实时信号测频方法以及锁相放大器
CN114371342B (zh) * 2022-03-21 2022-05-27 国仪量子(合肥)技术有限公司 Fpga及基于其的实时信号测频方法以及锁相放大器
CN118226123A (zh) * 2024-02-27 2024-06-21 青岛汉泰电子有限公司 一种用于示波器中快速查询fft峰值的方法

Also Published As

Publication number Publication date
CN106415310B (zh) 2018-02-02
EP3168638A4 (en) 2018-06-13
KR20170009878A (ko) 2017-01-25
EP3168638A1 (en) 2017-05-17
WO2016006079A1 (ja) 2016-01-14
US9857399B2 (en) 2018-01-02
KR101785714B1 (ko) 2017-11-06
JP5766896B1 (ja) 2015-08-19
JPWO2016006079A1 (ja) 2017-04-27
US20170205448A1 (en) 2017-07-20

Similar Documents

Publication Publication Date Title
US7555081B2 (en) Log-sampled filter system
Lim et al. A new algorithm for two-dimensional maximum entropy power spectrum estimation
CN106415310A (zh) 峰值频率检测装置、方法以及程序
AU751333B2 (en) Method and device for blind equalizing of transmission channel effects on a digital speech signal
JP3704336B2 (ja) 波形検出装置とそれを使用した状態監視システム
RU2751088C2 (ru) Способ и система обработки сейсмических данных
JP4274949B2 (ja) オーディオ・フィードバック処理システム
Cadzow Recursive digital filter synthesis via gradient based algorithms
Pavlenko et al. Identification of systems using Volterra model in time and frequency domain
CN107210046A (zh) 用于处理和分析信号的方法，以及实现这种方法的装置
Kay Maximum entropy spectral estimation using the analytical signal
CN107919134A (zh) 啸叫检测方法及装置和啸叫抑制方法及装置
Lim et al. A piloted adaptive notch filter
Zhou et al. A frequency band constrained filtered–x least mean square algorithm for feedback active control systems
CN107277669A (zh) 耳机的数字降噪滤波器生成方法及装置
Lim et al. FRM-based FIR filters with optimum finite word-length performance
JPH10503908A (ja) オーディオ信号の調性を決定するための方法および装置
Giampiccolo et al. A time-domain virtual bass enhancement circuital model for real-time music applications
CN109506762B (zh) 基于滤波器的水听器接收信号修正方法
CN109448747A (zh) 抗噪声波的生成方法、装置、计算机设备及存储介质
CN109788399A (zh) 一种音箱的回声消除方法及系统
CN108806711A (zh) 一种提取方法及装置
CN112116917B (zh) 基于相位跃变度的电抗器本体与风机声信号分离方法
CN102667508B (zh) 检测由检测器接收到的信号中的对应于事件的波前的方法
CN116755141B (zh) 一种深度域地震子波提取方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information
CB02 Change of applicant information

Address after: Chiba County, Japan

Applicant after: SFFT CO.,LTD.

Address before: Yokohama City, Kanagawa Prefecture, Japan

Applicant before: SFFT CO.,LTD.

GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20180202