JP5766896B1 - ピーク周波数検出装置、方法およびプログラム - Google Patents

ピーク周波数検出装置、方法およびプログラム Download PDF

Info

Publication number
JP5766896B1
JP5766896B1 JP2015509250A JP2015509250A JP5766896B1 JP 5766896 B1 JP5766896 B1 JP 5766896B1 JP 2015509250 A JP2015509250 A JP 2015509250A JP 2015509250 A JP2015509250 A JP 2015509250A JP 5766896 B1 JP5766896 B1 JP 5766896B1
Authority
JP
Japan
Prior art keywords
frequency
digital data
peak frequency
power
digital
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.)
Active
Application number
JP2015509250A
Other languages
English (en)
Other versions
JPWO2016006079A1 (ja
Inventor
一史 阿久津
一史 阿久津
Original Assignee
一史 阿久津
一史 阿久津
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 一史 阿久津, 一史 阿久津 filed Critical 一史 阿久津
Application granted granted Critical
Publication of JP5766896B1 publication Critical patent/JP5766896B1/ja
Publication of JPWO2016006079A1 publication Critical patent/JPWO2016006079A1/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R19/00Arrangements for measuring currents or voltages or for indicating presence or sign thereof
    • G01R19/04Measuring peak values or amplitude or envelope of ac or of pulses
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16Spectrum analysis; Fourier analysis
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16Spectrum analysis; Fourier analysis
    • G01R23/165Spectrum analysis; Fourier analysis using filters
    • G01R23/167Spectrum analysis; Fourier analysis using filters with digital filters
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16Spectrum analysis; Fourier analysis
    • G01R23/177Analysis of very low frequencies
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/02Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems using reflection of acoustic waves
    • G01S15/50Systems of measurement, based on relative movement of the target
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/35Details of non-pulse systems
    • G01S7/352Receivers
    • G01S7/356Receivers involving particularities of FFT processing
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/523Details of pulse systems
    • G01S7/526Receivers
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/14Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141Discrete Fourier transforms
    • GPHYSICS
    • G10MUSICAL INSTRUMENTS; ACOUSTICS
    • G10GREPRESENTATION OF MUSIC; RECORDING MUSIC IN NOTATION FORM; ACCESSORIES FOR MUSIC OR MUSICAL INSTRUMENTS NOT OTHERWISE PROVIDED FOR, e.g. SUPPORTS
    • G10G7/00Other auxiliary devices or accessories, e.g. conductors' batons or separate holders for resin or strings
    • G10G7/02Tuning forks or like devices
    • GPHYSICS
    • G10MUSICAL INSTRUMENTS; ACOUSTICS
    • G10LSPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/90Pitch 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)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Linguistics (AREA)
  • Human Computer Interaction (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Health & Medical Sciences (AREA)
  • Signal Processing (AREA)
  • Discrete Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (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

【課題】FFTにおける周波数解析の際に制約となる周波数分解能(f0)と時間窓長(T0)との間にある相反関係を回避し、所望の周波数分解能と所望の時間窓長でもって、信号波のピーク周波数の検出を可能にする。【解決手段】所定の周波数帯域(fcl〜fch)においてパワースペクトルが最大となるピーク周波数を検出するピーク周波数検出装置であって、デジタルデータ列の各要素をn乗(nは2以上の整数)するn乗部と、サンプリング周波数fsと周波数分解能ftgと時間窓長Ttgに応じて決まるN(Nは2のべき乗の整数)個のn乗されたサンプリング周波数fsのデジタルデータ列に対して高速フーリエ変換を行って得られるパワースペクトルの最大値に対応する周波数を仮想ピーク周波数として導出するFFT部と、前記仮想ピーク周波数を1/n倍した値をデジタルデータ列のピーク周波数として出力する1/n倍部と、を備え、n≧1/(ftg?Ttg)、fs/(n?ftg)≰N≰fs?Ttg、fs>2?n?fchを満たすピーク周波数検出装置。

Description

本発明は、パワースペクトルが最大となるピーク周波数の検出装置、方法およびプログラムに関する。
高速フーリエ変換(FFT)の時間窓長をT(s)、サンプリング周波数をf(Hz)、サンプリング数をN(Nは2のべき乗の整数)とすると、
Figure 0005766896
の関係がある。時間窓長の逆数f(Hz)は周波数分解能と呼ばれ、
Figure 0005766896
となる。
信号波をFFTで周波数解析をする場合、時間窓長はFFTを行う信号波の時間長さ、周波数分解能は信号波周波数の最小分解能、つまり周波数の検出精度となる。以下、この時間窓長のことを、「FFTにかける信号波の切り出し時間」「FFTにかける信号波の時間」などともいう。[数2]式の意味するところは、FFTの時間窓長と周波数分解能は相反関係にあるということである。この相反関係は、周波数分解能が小さい周波数のところで影響が顕著になる。例えば、0.01(Hz)の周波数分解能が必要であれば、その逆数の100(s)間の信号波をFFTにかけなければならない。また、FFTにかける信号波の時間が0.01(s)であれば、その逆数の100(Hz)の周波数分解能しか得られない。
この相反関係は、FFTを利用するにあたって制約となることがあり、相反関係の制約を回避すべく様々な工夫がなされてきた。例えば、受信した信号波をFFTにかけ、周波数領域に現れるスペクトルを構成するピーク周波数を求めるに際して、時間窓長Tを、予め決めた値に固定し、受信した信号波からT間隔で相互に重複しないようにずらしながら切り出した複数の期間のそれぞれについてFFTを適用し(即ち、複数回に亘ってFFTを適用し)、求めた複数のピーク周波数を平均化する、という方法が考えられる。この方法では、1回のFFTで求められる周波数分解能fよりも高い分解能(細かい周波数間隔)でピーク周波数を計算できることもある。しかし、FFTの操作回数を増やせば、FFTにかける信号波の切り出し時間がその回数分長くなると言う問題がある。
なお、デジタルデータ列にFFTを適用すると、FFTに供したサンプル数Nに対し、0からN−1までの一定の周波数間隔(周波数分解能(f)に等しい)での振幅の値が与えられるパワースペクトルが得られる。これらのパワースペクトルのうち、最大値のパワースペクトルに対応する周波数のことを、ピーク周波数と呼ぶ。最大値のパワースペクトルは、公知の方法により、例えば、1からN/2までのパワースペクトルを1個ずつ比較していけば特定できる。このとき、pの点でパワースペクトルが最大であれば、ピーク周波数fpkはfpk=p×fとなる。
特許文献1には、周波数分解能が12(Hz)以下、FFTにかける信号波の切り出し時間が10(ms)(水中の位置分解能7.5mに相当)以下を同時に満たしてピーク周波数を求める方法が説明されている(同文献の[0023]〜[0026]、[0089]〜[0090]参照。なお[0097]以降の数値例では、入力信号波からの切り出し時間が5(ms)に変わっているが、理由は不明である。)。周波数分解能12(Hz)を優先すると、[数2]式から、切り出し時間は、1/12(Hz)=83.3(ms)となる。一方、FFTを適用する場合にはそのサンプル数Nは2のべき乗であることが要請されるので、同文献では、信号切り出し時間を、N=1024に相当する102.4(ms)としている(同文献の[0055]参照)。FFTにかける入力信号波の切り出し時間を102.4(ms)とするためは、実際のデータが得られる時間5(ms)では不足するため、不足分の102.4−5=97.4(ms)間の部分にゼロ値データを付加し、このゼロ値データを加えて得られたデータ列(時間長102.4(ms)の部分)をFFTにかけて、ピーク周波数を求めている(同文献の[0103]参照)。
同文献記載の方法では、FFTにかけるデータ列の中に含まれる有効データ列は、全体の約1/10(同文献[0103]によれば、約1/20)に過ぎないので、入力信号波がアナログ信号波のときに一般に適用されるA/D変換をして得られるデジタル信号波にかけるデジタルバンドパスフィルタ(同文献図7では、デジタルBPF62)を、次数を高くかつ狭帯域にすると、デジタルバンドパスフィルタ通過後の信号が鈍ってさらに有効データ数が減るため、デジタルバンドパスフィルタのフィルタ特性を緩やかにしなければならない。つまり、同文献には、時間長の短い入力信号波の場合の外乱ノイズによるネガティブな影響を回避することが困難となるという潜在的な課題は何ら認識されておらず、従って、そのような潜在的な課題解決の手段は何ら示唆も教示もなされていない。
また、同文献の場合、FFTにかけるデジタル信号波(デジタルデータ列)の中には、FFTを適用する際に要請される時間窓長の約1/10(もしくは、約1/20)の長さの部分にしか実際のデータは存在しないために、多くのパワースペクトルが発生せざるを得ない。FFTを適用する実際の信号波に外乱ノイズが含まれている場合には、ピーク周波数の周辺にさらに余分なスペクトルが現れることになるので、ピーク周波数を特定することが困難になる可能性がある。
このように、入力信号のピーク周波数を求めるにあたり、相反関係の制約を回避するために種々の試みがなされているが、一般的な解決方法には未だ至っていない。
特開2012−247302号公報
本発明の目的は、FFTによる周波数解析の際に制約となる周波数分解能(f)と時間窓長(T)との間にあるf=1/Tという相反関係を回避し、所望の周波数分解能と所望の時間窓長でもって、信号波のピーク周波数の検出を可能にする装置、方法およびプログラムを提供することにある。
(1)上記目的を達成するためのピーク周波数検出装置は、所定の周波数帯域(fcl〜fch)においてパワースペクトルが最大となるピーク周波数を検出するピーク周波数検出装置であって、デジタルデータ列の各要素をn乗(nは2以上の整数)するn乗部と、サンプリング周波数fと周波数分解能ftgと時間窓長Ttgに応じて決まるN(Nは2のべき乗の整数)個のn乗されたサンプリング周波数fのデジタルデータ列に対して高速フーリエ変換を行って得られるパワースペクトルの最大値に対応する周波数を仮想ピーク周波数として導出するFFT部と、前記仮想ピーク周波数を1/n倍した値をデジタルデータ列のピーク周波数として出力する1/n倍部と、を備え、
n≧1/(ftg×Ttg
/(n×ftg)≦N≦f×Ttg
>2×n×fch
を満たす。なお本明細書において、fcl〜fchはfcl以上fch以下を意味するものとする。
周波数解析の対象となる信号波yが、単一の周波数f(Hz)を有する[数3]式で表される正弦関数であると仮定すると、yのn乗(nは2以上の正の整数)は[数4]式および[数5]式で表される。[数4]式においてnは奇数であって、[数5]においてnは偶数である。
Figure 0005766896
Figure 0005766896
Figure 0005766896
[数4]式および[数5]式は、yをn乗(n=2,3,4,・・・)すると、nが奇数か偶数かに関わらず、周波数成分n×f,(n−2)×f,(n−4)×f,…が表れるということを示している。
したがって、y=sin(2πft)をn乗して、そのn乗して得られた信号波から、周波数成分(n×f)の信号波を抽出し、その信号波のピーク周波数を求めることができれば、そこで求められた周波数を1/n倍すれば、それが、元の信号波の周波数fに対応することが理解できる。図1は、この関係を模式的に示したものである。
ここで、周波数解析に使用するFFTのサンプリング周波数をf、サンプル数をNとする。y=sin(2πft)をこのFFTにかけてピーク周波数を求めると、これがfの計算値であり、その時の周波数分解能fは、
Figure 0005766896
となる。
ここで、周波数成分が(n×f)の正弦関数y=sin(2πnft)について、yのピーク周波数を求めたときと同じ値のfとNを有するFFTにかけてピーク周波数を求めると、求めたyの周波数f'は、(n×f)の計算値となる。ここで注目すべきは、このyの周波数f'は、yについての周波数分解能と同じ周波数分解能fのもとで求められた、ということである。つまり、f'の周波数分解能もfである。したがって、f'を1/n倍したf'/nは、yの周波数fを、周波数分解能f/nで求めたものとなる。
したがって、yの周波数fを直接FFTにかけて求めたときの周波数分解能をf、時間窓長をTとし、yの周波数fをyから求めたときの周波数分解能をf、時間窓長をTとすると、
Figure 0005766896
Figure 0005766896
であり、
Figure 0005766896
の関係が成立する。
[数9]式について、n=2,3,4のときと、べき乗をしないときのf×T=1のグラフを図2に示す。
周波数分解能ftgは、元の信号波のピーク周波数を求めるにあたって、ユーザが希望するピーク周波数の周波数分解能である。時間窓長Ttgは、ユーザが希望するFFTの時間窓長である。ftgとTtgとは独立に定めることができる。
つまり、本発明は、元の信号波のピーク周波数を、以下の条件を満たして計算する方法である。
・ピーク周波数の周波数分解能f≦ftg
・FFTの時間窓長T≦Ttg
このためには[数9]式から、
Figure 0005766896
を満たす必要がある。
つまり、乗数nは、
Figure 0005766896
を満たさなければならない。
なお、nは2のべき乗、2,4,8,16,32,・・・に制限されず、2以上の整数から選べばよい。したがって、用途に応じて、必要で十分な最適な値を選択することができる。すなわち、本発明によると、FFTによる周波数解析の際に制約となる周波数分解能(f)と時間窓長(T)との間にあるf=1/Tという相反関係を回避し、所望の周波数分解能と所望の時間窓長でもって、信号波のピーク周波数を検出することが可能になる。具体的には、ftg×Ttg<1の範囲においても所望の周波数分解能と所望の時間窓長で信号波のピーク周波数を検出することができる。
信号波をn乗した後のデジタルデータ列をFFTにかける場合のピーク周波数の周波数分解能f、時間窓長Tには、[数9]式で示したとおりf×T=1/nの関係があり、これを図示すると図4のようになる。この図4に描かれたf×T=1/nのグラフにおける太線部分が、ftg、Ttgの要求を満たすf、Tとなる。
=ftgのときのTをTminとすると、Tの取りうる範囲は、
Figure 0005766896
となる。一方、
Figure 0005766896
だから、
Figure 0005766896
となる。
FFTの場合、Nは2のべき乗でなければならないため、Nの取りうる範囲は、
/(n×ftg)≦N≦f×Ttg(Nは2のべき乗)
となる。
もし、Nが存在しない場合は、Nが存在するまでnを大きくするか、fを大きくするか、もしくはその両方を行う。但し、FFTのサンプリング定理から、
Figure 0005766896
を満たしていなければならない。
該当するNが複数存在する場合、どのNを採用しても所望の周波数分解能ftgと所望の時間窓長Ttgを満たしている。但し、Nが小さい方がFFTの計算量は少なくなるため、通常は、
Figure 0005766896
に最も近いNを採用するのが良い。
このように、ftg×Ttg<1の範囲においても、所望の周波数分解能ftg、時間窓長Ttgでピーク周波数を検出できる乗数n、FFTのサンプリング周波数f、FFTのサンプル数Nを決めることができる。なお、乗数n、f、Nの決め方は、上記の手順によらなくても良い。適当に数値を当てはめて条件を満たすまで試行錯誤を繰り返したり、簡単なプログラムを作成して導出したり、どのような方法を用いても良い。
(2)上記目的を達成するためのピーク周波数検出装置は、前記所定の周波数帯域に含まれる周波数のデジタルデータ列を抽出する第一デジタルバンドパスフィルタを備え、前記n乗部には、前記第一デジタルバンドパスフィルタの出力が入力されてもよい。
この構成を採用すると、様々な周波数成分を含む信号から所定の周波数帯域の信号を抽出するために前段のアナログフィルタを簡略化または省略できるため、回路規模を小さくすることができる。
(3)上記目的を達成するためのピーク周波数検出装置は、サンプリング周波数fisのデジタルデータ列を1/r(rは2以上の整数)に間引いてサンプリング周波数をfにする間引き部を備え、前記第一デジタルバンドパスフィルタには前記間引き部の出力が入力されてもよい。
この構成を採用すると、ピーク周波数検出装置に入力されるデジタルデータ列のサンプリング周波数がfより大きい場合であっても、所望の周波数分解能ftg、時間窓長Ttgでピーク周波数を検出できる。
(4)上記目的を達成するためのピーク周波数検出装置は、受信デジタルデータ列をg倍(gは2以上の整数)に補間してサンプリング周波数をfにする補間部を備え、前記第一デジタルバンドパスフィルタには前記補間部の出力が入力されてもよい。
この構成を採用すると、ピーク周波数検出装置に入力されるデジタルデータ列のサンプリング周波数がfより小さい場合であっても、所望の周波数分解能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以上の整数)に間引いてサンプリング周波数をfとする間引き部を備え、前記FFT部には、前記間引き部の出力が入力されてもよい。
この構成を採用すると、ピーク周波数検出装置に入力されるデジタルデータ列のサンプリング周波数がfより大きい場合であっても、所望の周波数分解能ftg、時間窓長Ttgでピーク周波数を検出できる。
(7)上記目的を達成するためのピーク周波数検出装置において、前記特定の周波数帯域に含まれる周波数のデジタルデータ列を抽出する第一デジタルバンドパスフィルタと、前記n乗部の出力から第二の周波数帯域に含まれるデジタルデータ列を抽出する第二デジタルバンドパスフィルタとを備え、前記n乗部には、前記第一デジタルバンドパスフィルタの出力が入力され、前記FFT部には、前記第二デジタルバンドパスフィルタの出力が入力され、前記第二の周波数帯域は、ほぼn×fcl〜n×fchであってもよい。
(8)上記目的を達成するためのピーク周波数検出装置は、前記n乗部に代えて、入力されるデジタルデータ列をm乗するべき乗部(j)(j=1,2,…,k)と前記べき乗部(j)の出力から特定の周波数帯域fcl(j)〜fch(j)の信号を抽出するデジタルバンドパスフィルタ(j)とを備えるべき乗ブロック(j)をk(kは2以上の整数)段備える多重べき乗部を備え、
n=m×m×…×m
cl(j)≒(m×m×…×m)×fcl
ch(j)≒×m×…×m)×fch
であってもよい。
デジタルデータ列をべき乗する場合、乗数が小さいほど、べき乗することによって発生する余分な周波数成分を減らすことができる。余分な周波数成分が少ないほど、デジタルバンドパスフィルタで余分な周波数成分を除去しやすくなる。すなわち、1段のべき乗部と第二デジタルバンドパスフィルタで対応できる乗数nの上限よりも、べき乗部とデジタルバンドパスフィルタをユニット化して多段構成にして対応できる乗数nの上限の方が大きい。したがって、この構成を採用することにより、乗数として設定できるnの範囲を広くすることができる。
(9)上記目的を達成するためのピーク周波数検出装置は、ユーザの指示を受け付ける操作部と、前記指示に応じたn、f、Nの少なくともいずれかを設定するパラメータ設定部と、を備えてもよい。
この構成を採用することにより、次式を満たすn、f、Nをユーザが試行錯誤を繰り返して設定する必要がなくなる。
n≧1/(ftg×Ttg
/(n×ftg)≦N≦f×Ttg
>2×n×fch
なお請求項に記載された各手段の機能は、構成自体で機能が特定されるハードウェア資源、プログラムにより機能が特定されるハードウェア資源、又はそれらの組み合わせにより実現される。また、これら各手段の機能は、各々が物理的に互いに独立したハードウェア資源で実現されるものに限定されない。さらに、本発明は、方法としても、コンピュータプログラムとしても、コンピュータプログラムの記録媒体としても成立する。むろん、そのコンピュータプログラムの記録媒体は、磁気記録媒体であってもよいし光磁気記録媒体であってもよいし、今後開発されるいかなる記録媒体であってもよい。
本発明の実施形態にかかるスペクトルを示す図。 本発明の実施形態にかかるグラフ。 本発明の実施形態にかかるブロック図。 本発明の実施形態にかかるグラフ。 本発明の実施形態にかかるスペクトルを示す図。 本発明の実施形態にかかるブロック図。 本発明の実施形態にかかる波形図。 本発明の実施形態にかかるブロック図。 本発明の実施形態にかかるブロック図。 本発明の実施形態にかかるブロック図。 本発明の実施形態にかかる波形図。 本発明の実施形態にかかる波形図。 本発明の実施形態にかかるブロック図。 本発明の実施形態にかかるブロック図。 本発明の実施形態にかかるブロック図。 本発明の実施形態にかかる画面構成図。
以下、本発明の実施の形態を添付図面を参照しながら説明する。尚、各図において対応する構成要素には同一の符号が付され、重複する説明は省略される。
1.第一実施形態
第一実施形態では、サンプリング周波数fでサンプリングされた受信デジタルデータ列の、下限値fclと上限値fchとで定められた周波数帯域fcl〜fchにおけるピーク周波数fを、所望の周波数分解能ftgと所望の時間窓長Ttgを満たして検出するピーク周波数検出装置とこれを用いたドップラー計測器について説明する。
本発明の第一実施形態としてのピーク周波数検出装置1は、図3に示すように、第一デジタルバンドパスフィルタ(BPF)部11と、n乗部12と、第二デジタルバンドパスフィルタ(BPF)部13と、FFT部14と、1/n倍部15とを具備している。
はじめに乗数nと、周波数解析に使用するFFTのサンプリング周波数fと、周波数解析に使用するFFTのサンプリング数Nの決め方の手順を説明する。
(ステップ1.乗数nを定める)
まず、n≧1/(ftg×Ttg)の条件を満たすn(nは、2以上の整数)を定める。例えば、nを、n≧1/(ftg×Ttg)を満たす最小の整数とし、後述するステップ3でNが存在しなければ、nを1増加させて再計算をしても良いし、最初から大きめにしておいても良い。
(ステップ2.FFTのサンプリング周波数fを選択する)
FFTのサンプリング周波数fは、サンプリング定理を満たすため、
>2×n×fch
となるように選択する。なお、本実施形態では、ピーク周波数検出装置1に入力されるデジタルデータ列のサンプリング周波数がFFTのサンプリング周波数fとなる。fが回路上の制約等で上式を満たせない場合は、他の実施形態を使用する。
(ステップ3.FFTのサンプル数Nを定める)
次に、f/(n×ftg)≦N≦f×Ttg(Nは2のべき乗)を満たすNを選択する。もし、Nが存在しない場合は、Nが存在するまでnを大きくするか、fを大きくするか、もしくはその両方を行う。ただし、FFTのサンプリング定理から、
>2×n×fch
を満たしていなければならない。
以上説明した方法で選択される乗数n、f、Nは、n乗部12およびFFT部14に設定される。また、ピーク周波数を検出しようとする帯域の下限値fclと上限値fchは第一デジタルBPF11のカットオフ周波数として設定される。また第二デジタルBPF13のカットオフ周波数には、ピーク周波数を検出しようとする帯域の下限値fclと上限値fchのn倍の値が設定される。
(デジタルデータ列からピーク周波数を計算する)
乗数n、f、Nが設定されたピーク周波数検出装置1に対象とするデジタルデータ列を入力すると、以下のように、所望の周波数分解能ftgと所望の時間窓長Ttgを満たしてピーク周波数が検出される。
対象とするデジタルデータ列がピーク周波数検出装置1に入力されると、第一デジタルBPF11は、デジタルデータ列の幅広い周波数成分から、設定されている帯域外の余分な直流成分、低周波成分、高周波成分を排除し、単一の周波数fに近い周波数成分を抽出する。
ピーク周波数を検出しようとする帯域としてfcl〜fchが広すぎて、その間に大きなパワースペクトルが複数存在する場合、fcl〜fchを狭く設定し直して、図5のように、何回かに分けてピーク周波数を求めるのが良い。図5の例では、fcl〜fchが最初の周波数帯域、f、fが、大きなパワースペクトルを持つ2つの周波数、fcl'〜fch'が狭めた周波数帯域を示す。なお、fcl'〜fch'の間隔は一定である必要はなく、図5に示すように、周波数が高くなるほど、fcl'〜fch'の間隔を広げてもよい。
ここで第一デジタルBPF11を通過した後のデジタルデータ列を、A(1):a,a,a,…とする。A(1)には、周波数fのほかに、パワースペクトルの小さい若干の余分な周波数成分が含まれる。
次に、このA(1)をn乗部12に入力する。n乗部12では、A(1)の各要素を、上で定めたnを乗数としてべき乗する。n乗部12通過後のデジタルデータ列をB(n):b,b,b,…とすると、b=(a(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):c,c,c,…とする。C(n)には、周波数n×fのほかに、若干の余分な周波数成分が含まれる。
ここで、ピーク周波数検出装置1に入力されるデジタルデータ列のサンプリング周波数fを、f>2×n×fchとした理由は、C(n)の周波数帯域の上限がn×fchだからである。つまり、受信デジタルデータ列のサンプリング周波数fは、FFTのサンプリング定理を満たすことを条件とした。
次に、C(n)を、FFT部14に入力してピーク周波数を算出する。FFT部14では、C(n)に対し、設定されたサンプリング周波数fとサンプル数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が設けられており、変調器25には、局部発振回路26から局部発振周波数flocの信号が供給されている。変調器25の出力すなわち中間周波数信号は、アナログフィルタ27を介してアナログ/デジタル変換器(A/D)28に入力されてFFTに要請されるサンプリング周波数でデジタル信号に変換される。次にピーク周波数検出装置1に入力されてFFTにより周波数解析され、ピーク周波数検出装置1から、ピーク周波数が出力される。
本実施例では、以下の条件下で海水中の音波のピーク周波数を検出する方法を説明する。 海水中の音波伝搬速度C:1500m/s
送信周波数ftx:120kHz
局部発信周波数floc:137kHz
検出最大速度(水平方向)V:15m/s
検出速度精度(水平方向)V:0.15m/s
検出対象物の位置の精度(位置分解能)D:7.5m
また、音波は水平方向から斜め方向(θ=60°)に送受信するものとする。C≫Vの場合、ドップラー周波数fdopは、
Figure 0005766896
となる。検出最大速度Vが15m/sだから、
Figure 0005766896
となり、ドップラー信号は120±1.2kHzの範囲であり、観測される周波数帯域の幅Δfは、
Δf=2×1200Hz=2400Hz
となる。
ここで中間周波数fmidは、floc−ftx=137−120=17kHzとする。アナログフィルタ27は、
mid±(△f/2)=17000±1200Hz
の信号を通し、次段のA/D変換器28でエイリアシングが発生しないような特性とする。
検出速度精度Vが0.15m/sだから、
Figure 0005766896
となり、周波数分解能fは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のサンプリング周波数fを決定する。第一デジタルBPF11の周波数帯域fcl〜fchは、17000±1200Hzに合わせると、fcl=15800Hz、fch=18200Hzとなる。したがって、A/D変換器28のサンプリング周波数fは、
>2×n×fch=2×12×18200=436800Hz
となるため、f=510kHzとする。
FFTのサンプル数Nは、f/(n×ftg)≦N≦f×Ttgより、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は、
Figure 0005766896
となる。
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'17800=5.2Hzであり、±(ftg)/2=±6Hz以下に収まっている。なお、周波数分解能fのときのピーク周波数は、fの間隔で最も近い点となるため、周波数分解能fのときのピーク周波数の誤差は、±(f/2)以下となる。
本実施例において、第二デジタルBPF13から出力されるデジタルデータ列の実例を図7に示す。図7において、Tは、FFTにかける区間である。Tは、デジタルデータ列の振幅が落ち着くまでの区間である。Tの要素数は4096個であるので、Tの区間長は、4096/510000≒8.0msとなり、所望の時間窓長Ttg=10msを満たしている。
なお、区間Tで振幅が落ち着いた後のデジタルデータ列をFFTにかける方がより好ましい。Tの要素数は約600個であるので、Tの区間長は600/510000≒1.2msとなる。その場合でも、乗数nを大き目の12に設定したためピーク周波数を計算するために必要なデジタルデータ列の時間長は、1.2+8.0=9.2msであり、Ttg=10ms以下となっている。なお、本実施例では、T以降の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に入力されるデジタルデータ列のサンプリング周波数fが高く、乗数nと第一デジタルBPF11のカットオフ周波数fchがf>4×n×fchを満たす場合、第二実施形態を適用することが好ましい。
間引き部16では、デジタルデータ列を1/r(rは2以上の整数)に間引き、間引いた後のデジタルデータ列が、
>2×n×fch
となるようにする。fは間引き後のサンプリング周波数である。
1/rに間引くには、例えば、次のようにする。間引く前のデジタルデータ列を、
P(1):p,p,p,…とし、間引いた後のデジタルデータ列を、Q(1):q,q,q,…としたとき、q=p(r×i)(i=0,1,2,3,4,・・・)を行う。r=2のときは、q=p、q=p、q=p、q=pとなる。
間引く方法はこれ以外であっても良い。例えば、r=2のときに、q=(p+p)/2、q=(p+p)/2、q=(p+p)/2、…などとしても良い。
なお、r=2としたときに、
/(n×ftg)≦N≦f×Ttg (Nは2のべき乗)
を満たすNが存在せず、nを大きくすることもできない場合、本実施例は適用できない。この場合は第一実施形態を適用する。
以下、図6に示したドップラー計測器2に本実施例のピーク周波数検出装置3を適用した場合について説明する。なお、下記の数値データ以外の必要な数値データは、第一実施形態の数値例と同じとする。また、後の説明を簡略化するため、第一実施形態同様、f=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に間引くものとする。間引いた後のデジタルデータ列のサンプリング周波数をfとすると、
=fis/r=10.2MHz/20=510kHz>436.8kHz
となる。
間引く前のデジタルデータ列を、P(1):p,p,p,…とし、間引いた後のデジタルデータ列を、Q(1):q,q,q,…とすると、間引き部16では、q=p、q=p20、q=p40、q=p60、…を行うこととなる。
このQ(1)は、入力信号を、サンプリング周波数f=510kHzでA/D変換した信号列と同等であることが理解できる。従って、FFTのサンプリング周波数fは、f=510kHzとなる。
第一実施形態と同様、y=sin(2π17800t)をサンプリング周波数10.2MHzのA/D変換器28への入力信号と仮定し、擬似的なデータ列を作成した上で、上記の実施例の設定で実際に周波数解析を行なった結果は、第一実施形態の数値例と同じとなり、ftg、Ttgを満たす。
3.第三実施形態
図9は、本発明の第三実施形態を示すピーク周波数検出装置4の構成を示すブロック図である。ピーク周波数検出装置4は、第一実施形態の第二デジタルBPF13の後に間引き部17が追加された構成となっている。第一デジタルBPF部11、n乗部12、第二デジタルBPF13は、入力されるデジタルデータのサンプリング周波数が、第一実施形態ではfだったものが、fisとなっただけで、動作は同じである。受信デジタルデータ列のサンプリング周波数fが高く、乗数nと第一デジタルBPF11のカットオフ周波数fchに対して、f>4×n×fchの場合、計算量を削減するため、一般的には第二実施形態を適用するのが良い。しかし、本実施形態であっても、所望の時間窓長ftg、所望の時間窓長Ttgを満たした処理が可能である。
間引き部17では、受信デジタルデータ列を1/r(rは2以上の整数)に間引いて、間引いた後のデジタルデータ列のサンプリング周波数をfとしたとき、
>2×n×fch
となるようにする。間引いた後のデジタルデータ列は、FFT部14に入力され、その後の処置は、第一実施形態と同一である。
なお、r≧2でf/(n×ftg)≦N≦f×Ttg (Nは2のべき乗)を満たせない場合には本実施形態は適用できない。この場合は、第一実施形態を適用する。
以下、図6に示したドップラー計測器2に本実施例のピーク周波数検出装置4を適用した場合について説明する。なお、下記の数値データ以外の必要な数値データは、第一実施形態の数値例と同じとする。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に間引くものとする。間引いた後のデジタルデータ列のサンプリング周波数をfとすると、
=fis/r=10.2MHz/20=510kHz>436.8kHz
となる。FFTのサンプリング周波数fは、f=510kHzとなる。
第一実施形態の数値例と同様、y=sin(2π17800t)をサンプリング周波数10.2MHzのA/D変換器28への入力信号と仮定し、擬似的なデジタルデータ列作成した上で、上記の実施例の設定で実際に周波数解析を行なった結果は、第一実施形態の数値例と同じとなり、ftg、Ttgを満たす。
4.第四実施形態
10は、本発明の第四実施形態としてのピーク周波数検出装置5の構成を示すブロック図である。ピーク周波数検出装置5は、第一実施形態のデジタルBPF11の前に、補間部18が追加された構成となっている。補間部18を追加することにより、A/D変換後のデジタルデータ列を補間して増やすことによりサンプリング周波数を上げることができる。ピーク周波数検出装置5に入力されるデジタルデータ列のサンプリング周波数がfが低く、乗数nと第一デジタルBPF11の高いほうのカットオフ周波数fchに対して、2×fch<f<2×n×fchの場合、本実施形態を適用する。
補間部18では、受信デジタルデータ列をg倍(gは2以上の整数)に補間して、補間した後のデジタルデータ列のサンプリング周波数fが、f>2×n×fchとなるようにする。例えば、補間前のデジタルデータ列を、U(1):u,u,u,…とし、補間した後のデジタルデータ列を、V(1):v,v,v,…としたとき、補間部18は、
=u (i=0,1,2,3,・・・,(g−1))
=u (i=g,g+1,g+2,g+3,・・・,(2g−1))
=u (i=2g,2g+1,2g+2,2g+3,・・・,(3g−1))
・・・
となるように補間する。g=2のときは、v=u0、=u0、=u1、=u1、=u2、=u、…となる。補間する方法はこれ以外であっても良い。例えば、g=2のときに、
=u
=(u+u)/2
=u
=(u+u)/2
=u
=(u+u)/2
・・・
となるように補間しても良い。
以下、図6に示したドップラー計測器2に本実施例のピーク周波数検出装置5を適用した場合について説明する。なお、下記の数値データ以外の必要な数値データは、第一実施形態の数値例と同じとする。説明を簡略化するため第一実施形態同様、f=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倍に補間するものとする。補間した後のデジタルデータ列のサンプリング周波数をfとすると、
=fis×g=42.5kHz×12=510kHz>436.8kHz
となる。
補間する前のデジタルデータ列を、U(1):u,u,u,…とし、補間した後のデジタルデータ列を、V(1):v,v,v,…としたとき、補間部18は、
〜v11=u
12〜v23=u
24〜v35=u
36〜v47=u
・・・
となるように補間を行う。
補間後のデジタルデータ列V(1)の例を図で示すと、図11のように階段状になる。図11は、y=sin(2π17000t)をサンプリング周波数42.5kHzのA/D変換器28への入力信号と仮定し、擬似的なデジタルデータ列作成した上で、補間を行った例である。このV(1)のサンプリング周波数fは、510kHzなっている。従って、その後の処理は、図10の構成の場合と同様に、V(1)を、第一デジタルBPF11へ入力されるサンプリング周波数fの受信デジタルデータ列とみなして、同じ手法で周波数解析を行えばよい。
第一実施形態の数値例と同様、y=sin(2π17800t)をサンプリング周波数42.5kHzのA/D変換器28への入力信号と仮定し、擬似的なデジタルデータ列作成した上で、上記の実施例の設定で実際に周波数解析を行なった結果は、第一実施形態の数値例と同じとなり、ftg、Ttgを満たす。
ここで、n乗部12の前のデジタルフィルタを、帯域通過フィルタであるデジタルバンドパスフィルタにした理由を、図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に置き換えたものである。
多重べき乗部19は、図14に示すように、デジタルデータ列をm乗(mは2以上の整数)するべき乗部(j)(j=1,2,3,・・・,k)と、べき乗部(j)の出力から特定の周波数帯域fcl(j)〜fch(j)の信号を抽出するデジタルバンドパスフィルタ(j)(BPF(j))とから構成されているべき乗ブロック(j)がk段(kは2以上の整数)番号順に従属接続されたものである。
べき乗ブロック(j)のべき乗部(j)の乗数mは、置き換え前のn乗部の乗数nについて、
n=m×m×・・・×m
が成立するように選択される。
また、
cl(j)≒(m×m×…×m)×fcl
ch(j)≒(m×m×…×m)×fch
と設定する。
乗数mが偶数の場合、[数5]式より、yをm乗すると直流分が発生する。したがって、デジタルBPF(j)は、直流分を削除できる必要がある。
[数4]式、[数5]式で明らかなように、sin(2πft)をm乗する場合、乗数mを小さくすれば、発生する(m−2)×f以下の周波数成分の数を減らすことができる。そのため、その後のデジタルBPFで、必要な帯域以外の周波数成分を、より抑制しやすくなる。したがって、本実施形態を適用すると、より高い乗数nまで対応できるようになる。目安としては、n乗部の乗数nが16を越える場合は、本実施形態の採用を検討するのがよい。
本実施形態においても、第一実施形態同様、所望の周波数分解能ftgと所望の時間窓長Ttgを満たす。また本実施形態でも、近似計算を行わないため、計算するピーク周波数の精度劣化はない。
本実施形態では、べき乗部と第二デジタルBPF部を多段構成にしたため、第一実施形態と比べると計算量が増える。しかし、より高い乗数nまで、所望の周波数分解能ftgと所望の時間窓長Ttgを満たして精度劣化無くピーク周波数を計算できるという長所は、このマイナス面を補って余りあるものがある。
以下、図6に示したドップラー計測器2に本実施例のピーク周波数検出装置6を適用した場合について説明する。なお、下記の数値データ以外の必要な数値データは、第一実施形態の数値例と同じとする。後の説明を簡略化するため、第一実施形態同様、n=12、f=510kHz、N=4096、n=m×mとする。またm=4、m=3とする。
y=sin(2π17800t)をサンプリング周波数510kHzのA/D変換器28への入力信号と仮定して擬似的なデジタルデータ列作成した上で、第一デジタルBPF11に入力し、第一デジタルBPF11の出力を、A(1):a,a,a,…とする。
べき乗部(1)で、A(1)の各要素を4乗し、べき乗部(1)通過後のデジタルデータ列を、B(4):b,b,b,…
とする。すなわち、b= (a (i=0,1,2,3,4,…)となる。
このB(4)をデジタルBPF(1)に通し、デジタルBPF(1)デジタル通過後のデジタルデータ列をC(4):c,c,c,…とする。デジタルBPF(1)は、次数8次、カットオフ周波数63.2kHz(4×fcl)、72.8kHz(4×fch)のバタワース型のIIRフィルタとする。
次に、べき乗部(2)で、C(4)の各要素を3乗し、べき乗部(2)通過後のデジタルデータ列をD(12):d,d,d,…とする。すなわち、d =(ci(i=0,1,2,3,4,…)となる。
このD(12)をデジタルBPF(2)に通し、デジタルBPF(2)通過後のデジタルデータ列をE(12):e,e,e,…とする。デジタルBPF(2)は、次数8次、カットオフ周波数189.6kHz(4×3×fcl)、218.4kHz(4×3×fch)のバタワース型のIIRフィルタとする。
このE(12)について、サンプル周波数f=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/f=4096/510kHz=8.0ms<Ttg=10msであり、所望の時間窓長Ttgを満たしている。
6.第六実施形態
ここまで説明したピーク周波数検出装置1、3−6に、パラメータ設定部を追加してもよい。図15は、第一実施形態のピーク周波数検出装置1にパラメータ設定部20を追加した構成を示すブロック図である。パラメータ設定部20は、プロセッサ、メモリ、入出力機構を備えるコンピュータであって、図示しないキーボードやマウスやタッチパネルディスプレイ等の図示しない操作部を用いたユーザの入力に応じてピーク周波数検出装置1、3−6にパラメータの値を設定する。パラメータは、これまで説明したとおり、次のようなものである。
受信デジタルデータ列のサンプリング周波数fis
デジタルデータ列のサンプリング周波数f
所望の周波数分解能ftg
所望の時間窓長Ttg
n乗部の乗数n(nは2以上の整数)
第一デジタルBPFの周波数帯域 ほぼfcl〜fch(fcl<fch
第二デジタルBPFの周波数帯域 ほぼn×fcl〜n×fch
FFTのサンプリング周波数f
FFTのサンプリング数N
パラメータ設定部20は、上記全パラメータの数値をユーザが選択可能な入力値と対応付けて予めメモリに記憶しておき、ユーザの入力値に応じて設定してもよいし、一部のパラメータの数値をユーザが選択可能な入力値と対応付けて予めメモリに記憶しておき、ユーザの入力に応じて設定し、残りのパラメータは、入力された数値に応じてプロセッサが計算で求めても良い。
本実施形態では、周波数に対応した音階番号のみをユーザに入力させ、入力された音階番号に応じて全パラメータの数値をメモリから取得して設定する調律補助装置として使用されるピーク周波数検出装置について説明する。すなわち、パラメータ設定部20は、音階番号1〜88のいずれか一つの入力に対応するパラメータの数値をメモリから取得して設定する。図16は、音階番号を入力する状態におけるタッチパネルディスプレイの画面構成図の一例である。全パラメータの数値は、音階番号Pに応じて、前もって全て決められて不揮発性メモリに記憶されている。音階番号は平均律ピアノの鍵に対応し、音階番号1は最低音の周波数、音階番号88は最高音の周波数に対応する。すなわち、音階番号Pに対応する周波数fは次式[数21]によって定まる値が予め不揮発性メモリに記憶されている。
Figure 0005766896
そして例えば、P=49(f=442Hz)に対する各パラメータの数値は、以下のように定められている。
Figure 0005766896
なおftgは、fの2セント分にあたる周波数である。セント値は、周知の通り、2音間の周波数比を対数表現で表した値であって、100セントが平均律12音階の半音に相当する。Ttgは、時間窓長(1/ftg)の1/10である。fclは、fの50セント下の周波数である。fchは、fの50セント上の周波数である。Nは、f/(n×ftg)=24000/(16×0.5109)≒2936〜f×Ttg=24000×0.1957≒4697間の2のべき乗の整数でなければならない。
以上のように前もって取り決めておけば、パラメータ設定部20は、音階番号Pの入力に応じて全パラメータの数値を設定することができる。このように設定されたパラメータを用いることにより、受信デジタルデータ列のピーク周波数が計算できる。なお、ピーク周波数検出装置1に入力するデジタルデータ列は、図示しないマイク、A/D変換器を通じてピーク周波数検出装置1にリアルタイムで順次入力されても良いし、メモリに格納されていても良い。
P=49(f=442Hz)とし、y=sin(2π442t)をサンプリング周波数24kHzでサンプリングした擬似的なデジタルデータ列を作成した上で、実際に周波数解析を行なった結果は、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/f=4096/24kHz=0.171sであり、所望の時間窓長Ttg=0.195sを満たしている。なお、求めたピーク周波数の出力形式は、単位Hzでも良いし、fp=442Hzに対するセント値(+0.065セント)でも良い。
7.効果
以上説明した本発明の実施形態によると、FFTによる周波数解析の際に制約となる周波数分解能(f)と時間窓長(T)との間にあるf=1/Tと言う相反関係の問題を回避し、前記所望の周波数分解能ftgと前記所望の時間窓長Ttgでもって、信号波のピーク周波数の検出を可能にすることができる。受信デジタルデータ列のサンプリング周波数は、2×fch より大きければピーク周波数計算が可能となる。そして、近似計算やカーブフィッティング、平均化などの処理を必要としないため、ピーク周波数の計算精度の劣化は無い。
また、入力信号に余分な周波数帯域の成分(直流分、低周波分、高周波分)が含まれていても、S/N比が悪くても、殆ど振幅が無いような場合でも、第一デジタルBPF、第二デジタルBPFにより、ピーク周波数を計算するのに必要な帯域の周波数成分を抽出するため、問題なくピーク周波数を計算することができる。したがって、ピーク周波数検出装置の前段のハードウェアの要求仕様を下げることができ、小型化やコストダウンをはかることができる。具体的には例えば、第一デジタルBPFのカットオフ周波数fchに対し、入力信号のサンプリング周波数が2×fch より大きければピーク周波数が計算できるため、前段のアナログフィルタやA/D変換器などのハードウェアの要求仕様を下げることができる。また、べき乗部に設定する乗数nを大きめに設定すれば、デジタルバンドパスフィルタの次数を上げて急峻にし、外来ノイズに強くすることができる。
8.他の実施形態
尚、本発明の技術的範囲は、上述した実施例に限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々変更を加え得ることは勿論である。
例えば、IIR、バタワース型、8次のデジタルBPFを例示したが、他の形式のものも好適に使用できる。例えば、FIR(Finite impulse response)、チェビシェフ型などが好適に使用できる例として挙げられる。次数も8次に限らない。また、カットオフ周波数も、状況に応じて、抽出範囲を広くあるいは狭くすれば良い。また入力信号が、0を中心に振れる、高調波成分の殆ど無い信号(つまり、単一周波数のsinカーブに近い信号)であるなら、第一デジタルBPFをローパスフィルタにしても良い。
また反射エコーにおけるドップラー周波数の検出や調律補助装置への適用を例に挙げて本発明を説明したが、本発明の適用範囲はこれに限られるものではない。FFTにより信号波のピーク周波数の検出を行う用途に、本発明は広く一般的に適用可能である。
また、上記実施の形態に係る各機能部は、1または複数のLSI(Large Scale Integration)により実現されてもよく、また、複数の機能部が1のLSIにより実現されてもよい。また、集積化の手法としてはLSIに限定されるものではなく、専用回路または汎用プロセッサで実現することとしてもよい。さらには、LSI製造後にプログラムすることが可能なFPGA(Field Programmable Gate Array)や、LSI内部の回路セルの接続や設定を再構成可能なリコンフィギュラブル・プロセッサ(Reconfigurable Processor)を利用してもよい。また、半導体技術の進歩又は派生する別技術により、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は2以上の整数)するn乗部と、
    サンプリング周波数fと周波数分解能ftgと時間窓長Ttgに応じて決まるN(Nは2のべき乗の整数)個のn乗されたサンプリング周波数fのデジタルデータ列に対して高速フーリエ変換を行って得られるパワースペクトルの最大値に対応する周波数を仮想ピーク周波数として導出するFFT部と、
    前記仮想ピーク周波数を1/n倍した値をデジタルデータ列のピーク周波数として出力する1/n倍部と、
    を備え、
    n≧1/(ftg×Ttg
    /(n×ftg)≦N≦f×Ttg
    >2×n×fch
    を満たすピーク周波数検出装置。
  2. 前記所定の周波数帯域に含まれる周波数のデジタルデータ列を抽出する第一デジタルバンドパスフィルタを備え、
    前記n乗部には、前記第一デジタルバンドパスフィルタの出力が入力される、
    請求項1に記載のピーク周波数検出装置。
  3. サンプリング周波数fisのデジタルデータ列を1/r(rは2以上の整数)に間引いてサンプリング周波数をfにする間引き部を備え、
    前記第一デジタルバンドパスフィルタには前記間引き部の出力が入力される、
    ことを特徴とする請求項2に記載のピーク周波数検出装置。
  4. 前記デジタルデータ列をg倍(gは2以上の整数)に補間してサンプリング周波数をfにする補間部を備え、
    前記第一デジタルバンドパスフィルタには前記補間部の出力が入力される、
    ことを特徴とする請求項2記載のピーク周波数検出装置。
  5. n乗されたN個の前記デジタルデータ列から第二の周波数帯域に含まれるデジタルデータ列を抽出する第二デジタルバンドパスフィルタを備え、
    前記FFT部には、前記第二デジタルバンドパスフィルタによって抽出されたデジタルデータ列が入力され、
    前記第二の周波数帯域は、ほぼn×fcl〜n×fchである、
    請求項1に記載のピーク周波数検出装置。
  6. 前記第二デジタルバンドパスフィルタで抽出されたデジタルデータ列を1/r(rは2以上の整数)に間引いてサンプリング周波数をfとする間引き部を備え、
    前記FFT部には、前記間引き部の出力が入力される、
    請求項5に記載のピーク周波数検出装置。
  7. 前記特定の周波数帯域に含まれる周波数のデジタルデータ列を抽出する第一デジタルバンドパスフィルタと、
    前記n乗部の出力から第二の周波数帯域に含まれるデジタルデータ列を抽出する第二デジタルバンドパスフィルタとを備え、
    前記n乗部には、前記第一デジタルバンドパスフィルタの出力が入力され、
    前記FFT部には、前記第二デジタルバンドパスフィルタの出力が入力され、
    前記第二の周波数帯域は、ほぼn×fcl〜n×fchである、
    請求項1記載のピーク周波数検出装置。
  8. 前記n乗部に代えて、入力されるデジタルデータ列をm乗(mは2以上の整数)するべき乗部(j)(j=1,2,…,k)と前記べき乗部(j)の出力から特定の周波数帯域fcl(j)〜fch(j)の信号を抽出するデジタルバンドパスフィルタ(j)とを備えるべき乗ブロック(j)をk(kは2以上の整数)段備える多重べき乗部を備え、
    n=m×m×…×m
    cl(j)≒(m×m×…×m)×fcl
    ch(j)≒×m×…×m)×fch
    である請求項1〜4のいずれか一項に記載のピーク周波数検出装置。
  9. ユーザの指示を受け付ける操作部と、
    前記指示に応じたn、f、Nの少なくともいずれかを設定するパラメータ設定部と、
    を備える請求項1から8のいずれか一項に記載のピーク周波数検出装置。
  10. 所定の周波数帯域(fcl〜fch)においてパワースペクトルが最大となるピーク周波数を検出するピーク周波数検出方法であって、
    デジタルデータ列の各要素をn乗(nは2以上の整数)し、
    サンプリング周波数fと周波数分解能ftgと時間窓長Ttgに応じて決まるN(Nは2のべき乗の整数)個のn乗されたサンプリング周波数fのデジタルデータ列に対して高速フーリエ変換を行って得られるパワースペクトルの最大値に対応する周波数を仮想ピーク周波数として導出し、
    前記仮想ピーク周波数を1/n倍した値をデジタルデータ列のピーク周波数として出力する、
    ことを含み、
    n≧1/(ftg×Ttg
    /(n×ftg)≦N≦f×Ttg
    >2×n×fch
    を満たすピーク周波数検出方法。
  11. 所定の周波数帯域(fcl〜fch)においてパワースペクトルが最大となるピーク周波数を検出するピーク周波数検出プログラムであって、
    デジタルデータ列の各要素をn乗(nは2以上の整数)するn乗部と、
    サンプリング周波数fと周波数分解能ftgと時間窓長Ttgに応じて決まるN(Nは2のべき乗の整数)個のn乗されたサンプリング周波数fのデジタルデータ列に対して高速フーリエ変換を行って得られるパワースペクトルの最大値に対応する周波数を仮想ピーク周波数として導出するFFT部と、
    前記仮想ピーク周波数を1/n倍した値をデジタルデータ列のピーク周波数として出力する1/n倍部と、
    してコンピュータを機能させ、
    n≧1/(ftg×Ttg
    /(n×ftg)≦N≦f×Ttg
    >2×n×fch
    を満たすピーク周波数検出プログラム。
JP2015509250A 2014-07-10 2014-07-10 ピーク周波数検出装置、方法およびプログラム Active JP5766896B1 (ja)

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
JP5766896B1 true JP5766896B1 (ja) 2015-08-19
JPWO2016006079A1 JPWO2016006079A1 (ja) 2017-04-27

Family

ID=53888036

Family Applications (1)

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

Country Status (6)

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

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107402326B (zh) * 2017-07-20 2019-08-23 南京理工大学 一种改进s变换的有限窗长时频分析方法
CN107589424B (zh) * 2017-08-31 2021-12-17 努比亚技术有限公司 超声波采样方法、装置及计算机可读存储介质
CN114371342B (zh) * 2022-03-21 2022-05-27 国仪量子(合肥)技术有限公司 Fpga及基于其的实时信号测频方法以及锁相放大器

Family Cites Families (17)

* 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
KR970003966B1 (ko) * 1994-11-22 1997-03-24 삼성전자 주식회사 윈도우필터를 이용한 직접확산통신시스템의 수신기
JP3271504B2 (ja) * 1996-02-02 2002-04-02 三菱電機株式会社 周波数推定回路およびそれを用いたafc回路
US5799038A (en) * 1996-04-30 1998-08-25 Advantest Corporation Method for measuring modulation parameters of digital quadrature-modulated signal
JPH09298572A (ja) * 1996-04-30 1997-11-18 Advantest Corp Psk信号の搬送波周波数推定器
US7239674B2 (en) * 2003-06-04 2007-07-03 Honeywell Federal Manufacturing & Technologies, Llc Method of differential-phase/absolute-amplitude QAM
KR20070056818A (ko) * 2005-11-30 2007-06-04 주식회사 유컴테크놀러지 무선주파수인식 시스템
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 株式会社東芝 周波数同期回路、方法、プログラム及びこれらを用いた受信装置
JP5089460B2 (ja) * 2008-01-16 2012-12-05 三菱電機株式会社 伝搬遅延時間測定装置及びレーダ装置
US8390508B1 (en) 2010-04-05 2013-03-05 Raytheon Company Generating radar cross-section signatures
JP5094937B2 (ja) 2010-09-21 2012-12-12 三菱電機株式会社 周波数解析装置
JP2012247304A (ja) * 2011-05-27 2012-12-13 Sonic Corp 短時間信号のピークパワースペクトルを検出する方法及び装置
JP5840868B2 (ja) * 2011-05-27 2016-01-06 株式会社ソニック 周波数検出方法及び装置
JP5670836B2 (ja) 2011-05-27 2015-02-18 株式会社ソニック フーリエ変換でのサンプル数を削減した、短時間信号のピークパワースペクトルを検出する方法及び装置

Also Published As

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

Similar Documents

Publication Publication Date Title
EP2178082B1 (en) Cyclic signal processing method, cyclic signal conversion method, cyclic signal processing device, and cyclic signal analysis method
JP5670836B2 (ja) フーリエ変換でのサンプル数を削減した、短時間信号のピークパワースペクトルを検出する方法及び装置
AU2018208683B2 (en) Flow meter configuration and calibration
CN107210046B (zh) 用于处理和分析信号的方法,以及实现这种方法的装置
RU2011134419A (ru) Блок модулированных фильтров с малым запаздыванием
JP2013541024A (ja) 信号を解析し、瞬時周波数および短時間フーリエ変換を提供するための方法、ならびに信号を解析するためのデバイス
JP5766896B1 (ja) ピーク周波数検出装置、方法およびプログラム
JP5840868B2 (ja) 周波数検出方法及び装置
JP4734057B2 (ja) 移動物体検出装置
JP5492606B2 (ja) 演算装置、及び演算装置を備えた流量計
JP2012247304A (ja) 短時間信号のピークパワースペクトルを検出する方法及び装置
Huang et al. Resolution doubled co-prime spectral analyzers for removing spurious peaks
JP2000181472A (ja) 信号分析装置
JP6241131B2 (ja) 音響用フィルタ装置、音響用フィルタリング方法、およびプログラム
JP2002303645A (ja) 周波数計測装置、周波数計測方法およびレーダ装置
Jarrot et al. A time-frequency characterization framework for signals issued from underwater dispersive environments
White et al. Signal processing techniques
JP7056739B2 (ja) 波源方向推定装置、波源方向推定方法、およびプログラム
WO2021193637A1 (ja) 基本周波数推定装置、アクティブノイズコントロール装置、基本周波数の推定方法及び基本周波数の推定プログラム
KR100987306B1 (ko) 신호 처리 장치 및 그 방법
JP3964095B2 (ja) 大気温度測定方法および装置
JPH11160422A (ja) 電波高度計
JPH07231244A (ja) 狭帯域ディジタルフィルタ及びスペクトル分析器及びドップラ周波数推定器及び地殻同定器及び地磁気同定器及び地殻推定器及び地磁気推定器
JP5487062B2 (ja) 雑音除去装置
Nelson High-resolution correlation

Legal Events

Date Code Title Description
TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20150519

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20150617

R150 Certificate of patent or registration of utility model

Ref document number: 5766896

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

S111 Request for change of ownership or part of ownership

Free format text: JAPANESE INTERMEDIATE CODE: R313113

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350

S531 Written request for registration of change of domicile

Free format text: JAPANESE INTERMEDIATE CODE: R313531

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250