しかし、DSRVの場合、対象物である沈没潜水艦に接近するほど低速に制御され、ドップラソナーも高い速度分解能が要求されるにもかかわらず、短いパルスしか送信できないため速度分解能が低下するという問題点があった。たとえば、対象物との距離が1mの場合、反射波の到達遅れ時間が約1.3msであるため、サンプリング周波数が50kHzでもFFT点数は32点しかとることができず、周波数分解能は1.56kHzで検出速度分解能は約15knotとなり、超低速で接近する艦艇の速度を到底検出することができない。
そこで、送受波器を2つ設け、一方の送受波器で連続波を送信し、他方の送受波器で連続した反射波を受信することによって、長い信号から多くのサンプル点をとって、速度分解能を向上することが考えられる。しかし、単純に長い信号から長いサンプル点をとってFFT演算しようとすると極めて長時間の演算が必要になり、短時間の速度データの更新ができなくなるため、速度制御が不安定になるという問題点があった。
また、ドップラソナーはパルス信号の反射波遅延時間に基づいて対象物との距離を測定しているが、連続波ではこの距離測定ができないという問題点もあった。
この発明は、効率的な計算で超低速の速度分解能を得ることができる信号処理方法、信号処理装置、および、これらを備えたドップラソナーを提供することを目的とする。
この発明の信号処理方法は、所定の基準周波数付近に周波数スペクトルを有する信号を入力し、基準周波数をサンプリング周波数の4分の1の周波数にシフトするサンプリング周波数で信号をサンプリングし、サンプリングデータの周波数帯域を制限するFIRフィルタの各タップのフィルタ係数に対して信号の入力側から順に+1,−j、−1,+jまたは+1,+j、−1,−jの数列の任意の値から開始する繰り返しデータ列を乗算するとともに、各タップのフィルタ係数を4の倍数がデシメーションレートとなるように設定することで、FIRフィルタに入力される信号に対する周波数シフトと帯域制限とデシメーションとを同時に実行することを特徴としている。
また、この発明の信号処理方法は、所定の基準周波数付近に周波数スペクトルを有する信号を入力し、基準周波数をサンプリング周波数の4分の3の周波数にシフトするサンプリング周波数で信号をサンプリングし、サンプリングデータの周波数帯域を制限するFIRフィルタの各タップのフィルタ係数に対して信号の入力側から順に+1,+j、−1,−jまたは+1,−j、−1,+jの数列の任意の値から開始する繰り返しデータ列を乗算するとともに、各タップのフィルタ係数を4の倍数がデシメーションレートとなるように設定することで、FIRフィルタに入力される信号に対する周波数シフトと帯域制限とデシメーションとを同時に実行することを特徴としている。
また、この発明の信号処理装置は、所定の基準周波数付近に周波数スペクトルを有する信号を入力し、基準周波数をサンプリング周波数の4分の1の周波数にシフトするサンプリング周波数で信号をサンプリングするサンプリング手段と、サンプリングデータの周波数帯域を制限するように設定された各タップのフィルタ係数に対して、信号の入力側から順に+1,−j、−1,+jまたは+1,+j、−1,−jの数列の任意の値から開始する繰り返しデータ列を乗算するとともに、各タップのフィルタ係数を4の倍数がデシメーションレートとなるように設定することで、FIRフィルタに入力される信号に対する周波数シフトと帯域制限とデシメーションとを同時に実行するFIRフィルタと、を備えたことを特徴としている。
また、この発明の信号処理装置は、所定の基準周波数付近に周波数スペクトルを有する信号を入力し、基準周波数をサンプリング周波数の4分の3の周波数にシフトするサンプリング周波数で信号をサンプリングするサンプリング手段と、サンプリングデータの周波数帯域を制限するように設定された各タップのフィルタ係数に対して、信号の入力側から順に+1,+j、−1,−jまたは+1,−j、−1,+jの数列の任意の値から開始する繰り返しデータ列を乗算するとともに、各タップのフィルタ係数を4の倍数がデシメーションレートとなるように設定することで、FIRフィルタに入力される信号に対する周波数シフトと帯域制限とデシメーションとを同時に実行するFIRフィルタと、を備えたことを特徴としている。
この発明において、入力信号の周波数スペクトルは、注目領域の中心周波数である基準周波数付近に展開している。図1(A)は、この入力信号の周波数スペクトルの例を示す図である。この入力信号をサンプリング周波数fsでサンプリングするが、このサンプリング周波数fsを、図1(B)に示すように、サンプリング後の中心周波数fcがfs/4または3fs/4となるような周波数に設定する。このサンプリングには、たとえば、いわゆるアンダーサンプリングの手法を用いることができる。
このようにサンプリングすると、図1(C)に示すように、中心周波数がfcにバイアスされており、中心周波数を挟んで正負の範囲に周波数スペクトルが展開していてもその帯域幅がfc以下であれば、スペクトル同士がエリアシングで重なり合うことがない。また、このとき0Hz付近に生じる写像は1次写像とは限らず、2次写像などの多次写像である可能性があるが、中心周波数fcが事前に分かっているため、何次写像がどの付近に生じるかを予測することができ、最も利用しやすい0Hz付近の写像を用いることができる。なお、この技術はアンダーサンプリングのみならず、後述のデシメーションにおいても同様に適用することができる。
このように、サンプリング周波数fsでサンプリングされ、中心周波数fcの離散時間データ列となった信号は、
の指数関数列を乗算することによって中心周波数fcが0(DC)になるように周波数スペクトルをシフトすることができる。すなわち、データ数列x(n)のDFT変換から求まる周波数スペクトルが、
で求められるのに対し、データ数列x(n)に離散複素指数関数c(n)を乗算した周波数スペクトルXshift(k)が、
となることから周波数スペクトルX(k)が周波数軸に沿ってシフトされていることが分かる。すなわち、
によってスペクトルの注目領域の中心周波数fcを周波数ゼロとするように、スペクトル全体を周波数軸に沿ってシフトすることができる。
また、前記c(n)の指数部(−jΩc n)のnを、自然数Mを加算することによって(n+M)に置き換えた場合、すなわち,離散複素指数関数をM個シフトしてデータ数列に乗算した場合でも、
で明らかなように、周波数パワースペクトルはこのずれに影響されることなく同様にシフトされる。
そして、上述したようにサンプリング周波数fsと注目領域の中心周波数fcが、fc=fs/4またはfc=3fs/4となるような関係にサンプリングしていることにより、
Ωc =2π(fc/fs)=π/2
または、
Ωc =2π(fc/fs)=3π/2
となり、前記離散複素指数関数c(n)は、
となる。ここで、Ωc =π/2の場合を考えると、任意の整数値nに対して、
となり、+1,−j,−1,+jの4種類の値のみを取ることが分かる。
したがって、周波数スペクトルをシフトするために実際にx(n)とc(n)とを乗算する必要はなく、単にデータ数列x(n)を4個周期に、c(n)のnの値から簡単に割り出される+1,−j,−1,+jを乗算した場合に合わせて正負符号制御および実数虚数制御をするだけでよい。すなわち、c(n)がマイナス符号の場合には符号反転計算のみを行い、c(n)が実数の場合はx(n)の値を全て実数部として処理し、c(n)が虚数の場合はx(n)の値を全て虚数部として処理すればよい。なお、Ωc =3π/2の場合には、Ωc =π/2の場合と逆回りになり、+1,+j,−1,−jとなる。なお、+1→−j→−1→+jまたは+1→+j→−1→−jの繰り返しの先頭は+1,+j,−1,−jのうち任意のものでよい。
このように、fc=fs/4またはfc=3fs/4となるようなサンプリング周波数fsでサンプリングすることにより、サンプリングデータのサンプリング番号に基づいて符号制御および実数,虚数に割り振るのみの処理で周波数スペクトルのシフトを行うことができ、上記指数関数を実際に乗算して演算する必要がなくなるため、処理を大幅に簡略化することができる。
そして、この処理により周波数シフトされたサンプリングデータ列は、実数部のみのデータと虚数部のみのデータが交互に現れるため、後段のフィルタ演算などの演算においては、実数部の演算・虚数部の演算ともに通常の演算の1/2の演算量ですませることができる。すなわち、実数部の演算は、+1または−1が乗算されたサンプリングデータの実数部について行い、虚数部の演算は、+jまたは−jが乗算されたサンプリングデータの虚数部について行えばよく、図10や図13に示すように処理データ長の半分のデータ長の演算処理部でこれを実現できる。
また、上記のように中心周波数を0にシフトしたことにより、サンプリングデータの間引き(デシメーション)によって、図1(D)に示すように、0Hz(DC)を中心とした周波数スペクトルの引き延ばしが可能になる。サンプリングデータを間引きすることにより、長時間のサンプリングデータを少ないサンプル点(データ点数)で取り扱うことができる。
さらに、この発明では、デシメーションレートとして4n(n:正の整数)を用い、このデシメーションによって発生する折り返しスペクトルの重なり合い(エリアシング)を防止するためのローパスフィルタであるFIRフィルタのフィルタ係数に対して上記周波数シフトのための乗数である+1、−j、−1、+jまたは+1、+j、−1、−jを予め乗算している。デシメーションレートが4nであるから、入力されるサンプリングデータ列が4n進む毎にフィルタ演算が行われ、4個の周期で繰り返す+1、−j、−1、+jまたは+1、+j、−1、−jの乗数と同期し、同じサンプリングデータには常に同じ乗数が乗算されることになるため、FIRフィルタに入力するまえにサンプリングデータにこの乗数を乗算しておかなくても、フィルタ演算において同時に上記周波数スペクトル
のシフトを行うことができる。さらに、上述したように乗数として+1、−j、−1、+jまたは+1、+j、−1、−jを用いたことにより、サンプリングデータの実数部または虚数の一方が必ず0となるため、0となるタップ(係数演算)を省略することでフィルタ長を約1/2にすることができる。
これにより、周波数シフト、FIRフィルタ処理、およびデシメーションレート4nのデシメーション処理を一括して実行することができるとともに、フィルタ長を1/2にする(フィルタ演算量を1/2にする)ことができ、処理を大幅に簡略化することができる。
一方、有限長のサンプリングデータを切り出す場合、その周波数スペクトルを保存するためハニング窓などの窓関数を乗じてデータの切り出しを行うことがよく行われるが、この窓関数の各関数列に対して周波数スペクトルをシフトするための離散指数関数列c(n)を予め乗算しておき、これをサンプリングされたデータ列x(n)に乗ずることにより、窓関数によるサンプリングデータの切り出しと周波数スペクトルのシフトを同時に行うことができる。
さらに、この場合においても、図15に示すように上記周波数スペクトルのシフトと同様にサンプリング周波数fsと注目領域の中心周波数fcが、fc=fs/2またはfc=3fs/2となるような関係でサンプリングすることにより、窓関数列に対して+1,−j,−1,+jまたは+1,+j,−1,−jの乗数データ列を乗算したものをデータ列に乗算するのみで、簡略にサンプリングデータの切り出しと周波数スペクトルのシフトを同時に行うことができる。
この発明によれば、周波数スペクトルのシフト、周波数帯域の制限、デシメーションを一括した処理で行うことができるため、周波数スペクトルを有する信号の処理工程を簡略化することができる。特に、ディジタル化された周波数スペクトルのシフト、周波数帯域の制限、および、デシメーションを一括してFIRフィルタにおける処理で行うことができるため、信号処理を極めて簡略化することができる。
また、この発明によれば、ディジタル処理における周波数スペクトルのシフトを符号制御および実数部、虚数部への振り分けからなる簡略な処理で行うことができるため、周波数スペクトルを有する信号の処理をより簡略化することができる。
図2はこの発明の実施形態であるドップラソナーの概略構成図、図3は同ドップラソナーの受信部のブロック図である。このドップラソナーは、沈没潜水艦を救助する潜水艦救難艇(DSRV)に装着されるものであり、200ms毎に速度データを更新出力する。DSRVは、救助ハッチにドッキングする際は、沈没潜水艦に超低速で接近し、0.01ノットの精度で速度制御する必要がある。このため、このドップラソナーは0.01ノットの検出精度を有する。また、DSRVは救助を迅速に行うため、母船と沈没潜水艦との往復行程はある程度の高速で潜水・浮上する。このため、ドップラソナーも6ノット程度までの速度検出レンジを有している。
ドップラソナーは、上記検出精度および検出レンジを確保するため、沈没潜水艦とある程度以上の距離(20m〜25m以上)が離れているときはパルス信号を用いて速度計測を行い、沈没潜水艦と上記距離以下に接近したときは連続波信号を送信して反射波信号を長時間継続的に受信し、この長時間の信号を用いて速度計測を行う。上記パルス信号を用いて速度計測を行うモードをパルスドップラモードといい、連続波信号を用いて速度計測を行うモードを連続波ドップラモードという。
図2において、ドップラソナーは送信部1,受信部2および制御部3からなっている。送信部1は3方向に超音波ビームを送信する3つの送波器1aおよびこれら3つの送波器に対して高周波信号を印加する駆動回路などからなっている。駆動回路の動作は前記制御部3が制御する。
連続モードで動作する場合、前記送波器1aは継続的に超音波ビームを送信しているため、これを受波器として併用することができない。このため、受信部2は送波器1aとは別に、3個の受波器2aを備えている。各受波器2aは、前記3個の送波器1aにそれぞれ対応しており、対応する送波器と同じ方向に向けて設置されている。送波器1aが送信した超音波ビームは、沈没潜水艦などの対象物で反射して受波器2aに戻ってくる。この反射波信号は、このドップラソナーを搭載したDSRVと対象物との相対速度によるドップラ効果により周波数が遷移(ドップラシフト)する。このドップラシフト周波数を測定することにより、DSRVと対象物との相対速度を検出することができる。
また、パルスビームを送信した場合には、送信したのちその反射波を受信するまでの時間差を測定することができるため、この時間差に基づいてDSRVと対象物との距離を測定することができる。
前記受信部2は、3個の受波器2aのそれぞれについて、図3に示す受信回路を備えている。ただし、A/D変換器15およびDSP16と3系統で共有している。各送波器1a・受波器2aの組は、それぞれ3次元的に独立した方向に向けられているため、受信した反射波信号のドップラシフト周波数を測定し、このドップラシフト周波数から割り出された速度成分をベクトル合成することによってDSRVの前後,左右,上下方向の速度ベクトルを算出することができる。
なお、このドップラソナーが送信する超音波信号は400kHzであり、水中の音波の伝搬速度は1500m/秒であるため、超音波ビーム角を補正すると1ノット(=0.5144m/秒)に対応するドップラ周波数は、約100Hzである。
図3の受信回路において、受波器10(2aに相当)は対象物で反射した超音波信号を受信して電気信号に変換する。この反射波信号は、送信周波数とほぼ同じ約400kHzであるが、図4(A)に示すように、DSRVの移動速度に応じて0Hz〜±1kHz程度のドップラシフトを受けている。この反射波信号は、高周波アンプ11に入力される。高周波アンプ11はこの信号を後段の回路で処理可能なレベルまで増幅してミキサ12に入力する。ミキサ12には430kHzの変調信号が入力されており、この変調信号を前記反射波信号と混合する。この混合により、30kHz付近および830kHz付近に前記ドップラシフト周波数の写像スペクトルを生じる。バンドパスフィルタ13は30kHz付近に生じた下側うなりであるスペクトルのみを取り出してアンプ14に入力する。これにより400kHzのドップラシフト周波数信号が30kHzの中間周波にダウンコンバートされたことになる(図4(B)参照)。このダウンコンバートされたドップラシフト周波数信号は、下側うなり信号であるため受波器2aが受信したときの周波数スペクトルに対して周波数の高低が反転しているが、この反転は後段のDSP16における処理で再度反転させるか、または、反転しているものとして逆向きに処理することにより解消することができる。DSP16における反転処理は、上述の+1、−j、−1、+jまたはこれを反転した+1、+j、−1、−jのどちらで周波数スペクトルのシフトを行うかを選択することで行うことができる。
アンプ14はこの中間周波に変換されたドップラシフト周波数信号を適当なレベルまで増幅したのちA/D変換器15に入力する。A/D変換器15は、この信号を24kHzのサンプリング周波数でサンプリング(アンダーサンプリング)することによってディジタルデータに変換する。30kHzの信号を24kHzでアンダーサンプリングすることにより、30kHz±6kHzの周波数領域がサンプリングされ、この周波数領域のスペクトルが0(DC)〜12kHzおよび12kHz〜24kHzに現れる。したがって、中間周波において30kHz付近に生じていたドップラシフト周波数のスペクトルは、図4(C)に示すように、6kHz付近および18kHz付近に現れる。このようにドップラシフト周波数信号が、6kHzまたは18kHzにバイアスされたままサンプリングされることにより、ドップラシフトが送信周波数の下側に生じた場合でも実数値としてサンプリングすることができ、虚数側のA/D変換系統を設ける必要がなくなる。ここで、6kHz付近に現れるドップラシフト周波数のスペクトルは受信波信号に対してスペクトルが反転しており、18kHz付近に現れるドップラシフト周波数のスペクトルは受信波信号のスペクトルと同相である。
なお、このように30kHzの信号を24kHzでサンプリングすることにより、高い周波数の情報が失われるが、上述したようにDSRVの最大移動速度は6ノット程度であり、これによって生じるドップラシフトの周波数はDSRVの動揺による広がりを考慮しても1kHz程度であるため、上記±6kHzの周波数領域がサンプリングされることで十分DSRVの移動速度のレンジをカバーすることができる。
上記アンダーサンプリングによって、ディジタル変換されたドップラシフト周波数信号はDSP16に入力される。DSP16は、入力されたディジタル信号からドップラシフト周波数を算出し、そのドップラシフトを生じた超音波ビーム方向の速度成分を算出する。そして、前記3つの受波器2aの方向の速度成分をベクトル合成することによって、速度ベクトルを算出する。
DSP16は、最高精度の処理時においては、約1秒の時間幅のサンプリングデータを用いてFFTを行い、約1Hzの周波数分解能すなわち0.01ノットの速度分解能を得ているが、24kHzのサンプリング周波数で、そのまま1秒のサンプル点をとった場合には24×1024点=24576点という膨大なサンプル点でFFT演算する必要があり、演算量が膨大となって前記200ms毎のデータ更新期間内に処理することが困難になる。また、データ処理用のメモリを多く確保することが必要になり、大規模なDSPを使用する必要がある。一方、0.01ノットの速度分解能が必要な状態とは、対象物に対して極めて接近し速度を十分に落とした状態であるため、ドップラシフト周波数も十分低いと考えられる。このため、サンプリング周波数24kHzに対応する12kHzの検出可能最大周波数を確保する必要は全く必要なく、最大でも上記1kHz程度の検出可能最大周波数を有していれば、ドップラ周波数を検出可能である。したがって、このドップラソナーでは、サンプリングデータを間引き(デシメーション)してFFT処理をすることにより、検出可能最大周波数が低下するものの、少ないデータ数で演算量を増加させることなく周波数分解能を向上させている。
ここで、このDSP16におけるドップラシフト周波数の算出方法の概略を、図5,図6のフローチャートおよび図7〜図9を参照して説明する。なお、前記制御部3は、パルス信号の反射波を受信するまでの時間差により対象物(沈没潜水艦)までの距離を算出する機能を有している。
図5,図6の処理は200ms毎に繰り返し実行される速度算出処理動作を示している。この動作がスタートすると、まずs2で現在の動作モードを判断する。この動作モードは、パルスドップラモード(PWD),連続波ドップラモード(CWD)および距離測定モード(Temp.PWD)の3つであり、そのいずれかが前記制御部(メインCPU)3から指示される。制御部3は、図7に示す方式でパルスドップラモードと連続波ドップラモードを切り換える。すなわち、DSRVと対象物との距離が20m以下に接近したときパルスドップラモードから連続波ドップラモードに切り換え、両者の距離が25m以上に離れたとき連続波ドップラモードからパルスドップラモードに切り換える。このように、パルスドップラモードと連続波ドップラモードは、距離を基準にヒステリシスをもって
切り換えられる。
ここで、距離が接近したとき連続波ドップラモードに切り換えるのは、距離が接近していると送信信号を受信するまでの時間差が短いためパルスドップラモードでは長いパルスを送信することができず、FFT点数を多くとれないため速度分解能が低下するためであり、距離が離れているときパルスドップラモードに切り換えるのは、送受信の時間差が長いため長いパルスを送信でき、FFT点数を多くとり速度分解能を確保できるとともに、パルス送受信の時間差に基づいて並行して距離測定ができるためである。また、連続波ドップラモードでは超音波信号の送受信を並行して行うことから、送信信号の受信側への回り込みが問題となるが、距離が離れた対象物で反射した反射波信号は信号レベルが小さく上記回り込みによってカバーされるおそれがあるため、対象物との距離が離れた状態では連続波ドップラモードを用いることが困難であるためである。
上記のように、連続波ドップラモードでは、制御部3が距離測定をすることができないため、連続波ドップラモードが数十秒間継続するとパルスドップラモードの一種である距離測定モードに切り換えて距離を測定する。なお、このモード切換指示はこのDSP16のみならず送信部1にも与えられ、送信部1は、この指示に基づいて送波器1aから送信する超音波信号をパルス信号または連続波信号に切り換える。
対象物との距離が上記20m〜25m以上離れている場合には、制御部3からパルスドップラモードが指示される。この場合s3〜s10の処理を実行する。このパルスドップラモードの処理動作は一般的なドップラソナーの動作とほぼ同じである。s3では速度計測をするための前処理を実行する。前処理としては、A/D変換器15から入力されるサンプリングデータのレベルに基づいて反射波信号がどこに存在するかを検索し、FFT用のデータとしてそのパルス幅に応じた32〜1024点のデータを切り出す処理や後述の連続波ドップラモードで使用される処理サイクル数カウンタcycle#を1にリセットする動作などが含まれる。上記切り出されたサンプリングデータを高速フーリエ解析(FFT)して周波数スペクトルに展開し(s4)、この周波数スペクトルのピークを検索することによって、ドップラシフト周波数を割り出す(s5)。こののち過去に割り出されたドップラシフト周波数との平均化などの後処理を行って(s6)、今回の処理サイクルにおけるドップラシフト周波数を確定する。そして、このドップラシフト周波数が適正な範囲の値であるかをs7で判断する。すなわち、このDSRVは最大6ノット程度の速度でしか移動しないものであるため、ドップラシフト周波数は船体の動揺などによるマージンを考慮したとしても±1kHz程度である。したがって、±1kHzをしきい値とし、割り出されたドップラシフト周波数がこのしきい値を超える場合は計測が正常に行われていないとしてエラー信号を出力する(s10)。割り出されたドップラシフト周波数が正常な範囲の値であれば、上述した3系統の速度成分をベクトル合成することによってこのDSRVの速度ベクトルを算出し(s8)、これを他の航行制御装置に対して出力する(s9)。
一方、s2において、制御部3から指示された動作モードが連続波ドップラモードである場合には、s12以下の動作に進む。s12ではこの連続波ドップラモードの連続処理サイクル数をカウントする処理サイクル数カウンタcycle#の内容をチェックする。cycle#=1、すなわち、今回が連続波ドップラモードに切り換わった最初の処理である場合にはs12からs13に進む。s13では、現在蓄積されているデータ数に基づいてデシメーションするか否かを判断する。すなわち、連続波ドップラモードへの切り換えは、制御部3がこのDSP16および送信部1に指示するものであり(実際には制御部3が共有メモリにモード指示を書き込み、送信部1およびDSP16はそれを参照しにゆくことで指示を受け取る)、200ms毎に実行されるこの処理においてどの程度のデータが蓄積されているかはA/D変換されたサンプリングデータのメモリアドレス用のカウンタ値を参照することによってどこからが連続波ドップラモードであるかが判断される。一般的に連続波ドップラモードへの直後にFFTに使用可能なサンプリングデータ数は図8に示すように増加する。ここで、切換直後一定時間(100ms程度)のデータは無効なデータとしてオミットするようにしている。連続波ドップラモードに切り換えられた第1回目の処理においては、デシメーションするほどのデータが蓄積されていない場合が多いため、このような場合にはs13からs14に進んで蓄積されているデータのうち最新の1024点を切り出し、このデータをそのまま用いてFFTを行う(s19)。一方、1回目の処理が行われるとき、既に2048を超える有効なサンプリングデータが蓄積されている場合には、デシメーションレート=4(4個のデータ毎に1つのデータを用いて3個のデータ間を引くこと、以下同じ)、FFT点数=512と決定する(s15)。そしてこの処理に必要な2048個のデータ列を切り出してデモジュレーションにより周波数スペクトルをシフトする(s16)。このデモジュレーションは、上述したように、切り出されたサンプリングデータに離散複素指数関数列を乗算することによって行われる。
この処理で周波数シフトされたサンプリングデータ列に対してローパスフィルタ処理を実行する(s17)。これは、こののちデシメーションによって実質的にサンプリング周波数を低下させるため、高い周波数帯域にスペクトルがあるとエリアシングによりスペクトルが乱れることを防止するためである。そして、この周波数シフトされたデータ列からデシメーションレートのデータ数毎にFFT点数のサンプリングデータを抽出し(s18)、このデータを用いてFFTを実行する(s19)。
FFTによりサンプル数の周波数分解能で周波数スペクトルに分解したのち、そのピークを検索する(s20)。検出されたピーク周波数が上記しきい値(1kHz)を超える不適当なものであるかを判断し(s21)、不適当な値であればエラー信号を出力する(s26)。適当な範囲の値であれば過去の測定値の影響を残すために、平均化やIIRフィルタリングなどの処理を行って(s22)、計測精度を安定化させる。たとえば、パルスドップラモードから連続波ドップラモードに移行したとき、両モードで割り出されたピーク周波数にギャップが存在した場合などに、他の航行制御装置の制御が不安定になるのを防止するためこの処理が行われる。図9に示すcycle#≦6までの移行期には、パルスドップラモード時に検出されたピーク周波数の影響を残して出力する速度ベクトルを安定化する。このためにIIRフィルタを用いる場合には、例えば、
y(n)=(1−α)x(n)+αy(n−1)
の特性のものを用いればよい。ここで、y(0)を連続波ドップラモードに移行する直前のパルスドップラモードでのドップラシフト周波数、x(n)を連続波ドップラモードでのドップラシフト周波数とする。このフィルタ関数H(z)は、
H(z)= (1−α) / (1−αz-1)
で表される。連続波ドップラモードに移行した直後の、cycle#=1〜6までの約1秒間の処理サイクルにおいてパルスドップラモードにおける測定値の影響を残すため、計数αは0.5〜0.7の範囲で選択する。α=0.5のときcylce#=6のときのy(0)の影響は1.56%、α=0.6のとき5.00パーセント、α=0.7のとき約11.76%となる。また、逆に連続波ドップラモードからパルスドップラモードに移行する場合にも同様の処理を行えばよいが、αの設定値は上記の範囲に限定されるものではなく、他の値の範囲をとりうる可能性がある。
なお、ここではパルスドップラモードの測定値を最後の1個のみIIRフィルタ処理にかけているが、どの距離の測定時においても複数回の測定値を平均化して出力される速度データを安定化するようにしてもよい。
上記平均化・フィルタリング処理ののち3組の送波器1a,受波器2aによる3系統の計測結果をベクトル合成して速度ベクトルを算出する(s23)。算出された速度ベクトルを他の航行制御装置に対して出力する(s24)。以上の処理ののち処理サイクル数カウンタcycle#に1を加算して(s25)、リターンする。
DSRVが対象物に接近した状態であると、制御部3から継続的に連続波ドップラモードが指示され、s2→s12以下の処理が繰り返し実行される。cycle#=2の処理サイクルにおいては、s2→s12→s28→s29に進む。図8,図9に示すように、連続波ドップラモードに移行したのち処理サイクルを繰り返すと、蓄積される連続波のサンプリングデータ数も増加するため、cylc#の値に応じてFFTに使用するサンプリングデータ点数を切り換える。図8に示すようにcycle#=2のタイミングにおいては7200程度のデータ点数が使用可能になっているため、デシメーションレートを4に設定し、FFT点数を1024に設定する(s29)。これにより、4096サンプルの時間幅すなわち約170msの時間幅のデータに基づいて速度を割り出すことができ、約5.86Hzの周波数分解能を得ることができる。これに基づいて速度ベクトルを算出すれば、約0.06ノットの速度分解能を得ることができる。s29におけるデシメーションレートおよびFFT点数の設定ののちs16に進む。
また、cycle#=3の処理サイクルにおいては、s2→s12→s28→s30→s31に進む。図8に示すようにcycle#=3においては約12000の有効データが蓄積記憶されているため、デシメーションレート=8、FFT点数1024に設定する(s31)。これにより、8192サンプルの時間幅すなわち約340msの時間幅のデータに基づいて速度を割り出すことができ、約2.93Hzの周波数分解能を得ることができる。これに基づいて速度ベクトルを算出すれば、約0.03ノットの速度分解能を得ることができる。s31におけるデシメーションレートおよびFFT点数の設定ののちs16に進む。
さらに、cycle#=4の処理サイクルにおいては、s2→s12→s28→s30→s32に進む。s32では前の処理サイクルの測定結果に基づいて今回のピーク周波数を推定する。このピーク周波数の推定はDSRVが200msの短時間に急激に加速・減速できないこともに基づき、前回測定された速度から加速・減速できる範囲のドップラシフト周波数が全てこの350Hz内に含まれるかで判断される。この推定ののちs33からs34に進み、この推定されたピーク周波数が350Hz以下であるか否かを判断し、ピーク周波数が350Hz以下であると推定される場合にはs35に進み、ピーク周波数が350Hzを超えると判断された場合にはs36に進む。s35では、フィルタのカットオフ周波数を350Hzに設定するとともにデシメーションレート=16、FFT点数
1024に設定してs16に進む。上記フィルタのカットオフ周波数はs17のフィルタリング処理において適用される。デシメーションレートを16に設定すると24kHzでサンプリングされたデータが実質的に1.5kHzのサンプリング周波数に低下するためエリアシングを防止する必要からこのようなカットオフ周波数に設定する。デシメーションレート=16、FFT点数1024に設定することにより16384サンプルの時間幅すなわち約670msの時間幅のデータに基づいて速度を割り出すことができ、約1.46Hzの周波数分解能を得ることができる。これに基づいて速度ベクトルを算出すれば、約0.015ノットの速度分解能を得ることができる。
一方、s34で今回のピーク周波数が350Hzを超える、すなわち、DSRVが3.5ノットを超える速度で移動していると判断された場合にはs36に進んで、フィルタのカットオフ周波数を800Hzに設定するとともにデシメーションレート=12、FFT点数1024に設定してs16に進む。上記フィルタのカットオフ周波数はs17のフィルタリング処理において適用される。デシメーションレートを12に設定すると24kHzでサンプリングされたデータが実質的に2kHzのサンプリング周波数に低下するためエリアシングを防止するためにこのようなカットオフ周波数に設定する。デシメーションレート=12、FFT点数1024に設定することにより12288サンプルの時間幅すなわち約500msの時間幅のデータに基づいて速度を割り出すことができ、約2Hzの周波数分解能を得ることができる。これに基づいて速度ベクトルを算出すれば、約0.02ノットの速度分解能を得ることができる。
また、cycle#=5の処理サイクルにおいては、上記cycle#=4のときと同様にs2→s12→s28→s30→s32に進み、前の処理サイクルの測定結果に基づいて今回のピーク周波数を推定する。そしてs33→s37→s38に進み、この推定されたピーク周波数が350Hz以下であるか否かを判断する。ピーク周波数が350Hz以下であると推定される場合にはs39に進み、ピーク周波数が350Hzを超えると判断された場合には前記cycle#=4のときと同様にs36に進む。s39では、フィルタのカットオフ周波数を350Hzに設定するとともにデシメーションレート=20、FFT点数1024に設定してs16に進む。デシメーションレート=20、FFT点数1024に設定することにより20480サンプルの時間幅すなわち約850msの時間幅のデータに基づいて速度を割り出すことができ、約1.2Hzの周波数分解能を得ることができる。これに基づいて速度ベクトルを算出すれば、約0.012ノットの速度分解能を得ることができる。
そして、図8および図9に示すようにピーク周波数が350Hz以下であれば、cycle#=6以後の処理サイクルにおいては、24000点以上のサンプリングデータ点数を使用することができるため、最高分解能でピーク周波数を割り出す。まず、上記cycle#=4のときと同様にs2→s12→s28→s30→s32に進み、前回の処理サイクルの測定結果に基づいて今回のピーク周波数を推定する。そしてs33→s37→s40に進み、この推定されたピーク周波数が350Hz以下であるか否かを判断する。ピーク周波数が350Hz以下であると推定される場合にはs41に進み、ピーク周波数が350Hzを超えると判断された場合にはs42に進む。s41では、フィルタのカットオフ周波数を350Hzに設定するとともにデシメーションレート=24、FFT点数1024に設定したのちs43に進む。デシメーションレート=24、FFT点数1024に設定することにより24576サンプルの時間幅すなわち約1秒の時間幅のデータに基づいて速度を割り出すことができ、約1Hzの周波数分解能を得ることができる。これに基づいて速度ベクトルを算出すれば、目標である約0.01ノットの速度分解能を得ることができる。
一方、s34で今回のピーク周波数が350Hzを超える、すなわち、DSRVが3.5ノットを超える速度で移動していると判断された場合にはs42に進んで、フィルタのカットオフ周波数を800Hzに設定するとともにデシメーションレート=12、FFT点数1024に設定したのちs43に進む。デシメーションレート=12、FFT点数1024に設定することにより12288サンプルの時間幅すなわち約500msの時間幅のデータに基づいて速度を割り出すことができ、約2Hzの周波数分解能を得ることができる。これに基づいて速度ベクトルを算出すれば、約0.02ノットの速度分解能を得ることができる。
なお、処理サイクルを繰り返すことによってcycle#の値が限りなく増加することを防止するためs41,s42においてcycle#の値を6に固定している。
s43では、現在トラッキングモードであるかを判断する。トラッキングモードとは、これまで測定されたスペクトル形状からドップラシフト周波数信号の帯域幅が狭い範囲に限定されており、且つ、継続的に測定値が安定しており、今回も測定値に殆ど変化がないと考えられる場合のモードである。このような場合には、デモジュレーションを省略して直接デシメーションを行うことができる。たとえば、24kHzのサンプリングデータをデシメーションレート24でデシメーションすると実質的に1kHzのサンプリングデータとなり、6kHz付近に展開しているドップラシフト周波数のスペクトルは、各所に折り返し写像スペクトル(エイリアス)を生じる。しかし、ドップラシフト周波数のスペクトルが狭い帯域幅に限定されており、他のスペクトルがない場合には、これらのエイリアスが重なることはなく、いずれかのエイリアスをデシメーションされた周波数スペクトルとして取り扱ってFFTをすることが可能になる。24kHzのサンプリングデータをデシメーションレート24または12でデシメーションした場合、0Hz(DC)を中心周波数とするエイリアスも生じるため、これを処理用の周波数スペクトルとして用いることにより、通常処理と同様のFFTが可能になる。
s43で、現在トラッキングモードであると判断された場合にはs44に進み、周波数スペクトルの帯域幅が十分に狭く限定されているかを判断する。帯域幅が十分に狭く限定されていると判断された場合には、そのままs18に進む。周波数スペクトルの帯域幅がデシメーションによって重なり合う程度に広い場合には、ピーク周波数を求めるために不要な周波数帯域をカットするためバンドパスフィルタの処理を行ったのち(s45)、s18に進む。
s43でトラッキングモードでない場合と判断された場合には、s16に進み、通常どおりの処理を実行する。
また、連続波ドップラモードの継続中に20秒に1回の頻度で、制御部3から距離測定モードが指示される。この場合、DSP16は何の処理も行わずにそのままリターンする。このモードにおいて、送信部1は連続波ドップラモードの継続中に連続超音波信号の送信を一瞬停止して短時間のパルスを送信し、制御部3が反射波の時間遅れに基づいて対象物との距離を測定する。この距離測定処理は前記制御部3が実行するため、このDSP16は処理を行う必要がなく、また、このとき連続波ドップラモードの速度測定もできないため、この間DSP16はフリーズする。この距離測定モードによって測定された距離が図7のモード切換しきい値よりも短い場合には連続波ドップラモードに復帰し、距離が範囲外であれば通常のパルスドップラモードに切り換える。
連続波ドップラモードに復帰したとき、連続波ドップラモードを再開した時点から、新たにcycle#=1から処理サイクルを再開するようにしてもよく、バッファ内にストアされている以前のA/D変換データと再開後に収集されたデータとを連続したデータと見なして測定を行うようにしてもよい。この場合には、中断していた時間を考慮して中断前のデータと再開後のデータを滑らかに接続する。接続の方法としては、補間関数を用いて補間する方法、前後のデータから自己相関に基づいて接続する方法などがある。
さらに、パルス波による距離測定モード時においてもカウンタをフリーズすることなくデータ収集を継続するようにしてもよく、また、カウンタをフリーズさせずにメモリのアドレスをカウントアップしながら0などのブランクデータを書き込むようにしてもよい。この場合でも、パルス波の距離測定モードの間のデータは所定の方法で補間すればよく、また、若干精度がラフになるものの補間することなく中断区間のあるデータをそのまま測定に用いてもよい。また、大型船に用いられるドップラソナーの場合には、パルスモード時に収集したデータでドップラシフト(速度)を求めることもできる。
以上のような動作により、対象物と距離が離れている場合には、パルス信号を用いて速度計測を行い、対象物との距離が近い場合には、連続波信号を用いて速度測定を行うことにより、必要な速度分解能を確保しつつ、パルス信号のみを用いた場合の問題点である近距離における速度分解能の低下を解決するとともに、連続波のみを用いた場合に同時に速度と距離を測定することができないなどの問題点を解決した。
ここで、上記フローチャートでは、DSP16は、デモジュレーション、フィルタリング、デシメーションを順次行ってFFT用のデータを準備するようになっているが、実際には、上記手順を順次行うのではなく、デモジュレーション・フィルタリング・デシメーションを一気に行う簡略化された演算処理によってこれを行い、さらに、過去の処理済データを利用することによってより処理時間を短縮している。以下、この簡略化された演算処理方式について説明する。
サンプリング周波数fsによってA/D変換されたデータ数列x(n)から得られる周波数スペクトルの注目領域(スペクトルが展開している周波数帯)を拡大するためにデシメーション処理を行う。
このデシメーション処理を可能にするために前記注目領域の中心周波数(6kHz)をゼロ周波数(DC)にシフトする。すなわち、周波数軸に対してスペクトルを並行移動する。この操作は、上述したように、前記サンプリングデータ数列x(n)の各データに対して離散複数指数関数列c(n)を乗算することによって行う。
サンプリング周波数24kHzでサンプリングされ、中心周波数6kHzの離散時間データ列となった信号に対して、〔数8〕の指数関数列を乗算することによって中心周波数6kHzが0(DC)になるように周波数スペクトルをシフトする。この離散複素指数関数列は、複素単位乗数データ列(+1、−j、−1、+jまたは+1、+j、−1、−jの任意の値から開始する数列)の4種類の値のみを取るから、実際の処理では乗算を行う必要がなく、c(n)がマイナス符号の場合には符号反転計算のみを行い、c(n)が実数の場合はx(n)の値を全て実数部とし、c(n)が虚数の場合はx(n)の値を全て虚数部とするのみである。
図10は、A/Dデータバッファx(n)から複素数バッファX(n)への転記方式を説明する図である。A/Dデータバッファx(n)は、アンダーサンプリングされたサンプリングデータ(A/Dデータ)列を記憶するバッファであり、複素数バッファX(n)は、中心周波数=0に周波数スペクトルシフトされた虚数部を含むサンプリングデータ列を記憶するバッファである。
この図において、x(0)はそのままX(0)の実数部に転記され、X(0)の虚数部には0が書き込まれている。x(1)は正負の符号を反転されたのちX(1)の虚数部に転記され、X(1)の実数部には0が書き込まれている。x(2)は正負の符号を反転されたのちX(2)の実数部に転記され、X(2)の虚数部には0が書き込まれている。x(3)はそのままX(3)の虚数部に転記され、X(3)の実数部には0が書き込まれている。このように、離散複素指数関数の演算結果を複素単位乗数データ列(+1、−j、−1、+jまたは+1、+j、−1、−jの任意の値から開始する数列)の値にしたがって順次符号反転および転記を繰り返すのみでこの周波数シフトを行うことができ、指数関数を実際に乗算して演算する必要がなくなり、処理を大幅に簡略化することができる。
さらに、前記複素数バッファX(n)の実数部Real(n)、Imaginary(n)のうち一方は必ず0であるため、上記規則に基づいて0になる側が分かっていれば0を記憶するバッファを省略してバッファの記憶領域を実質的に半分にすることができる。
このように処理が簡略化されたとはいえ、連続波ドップラモード時にFFT処理のたびにデシメーション前のA/D変換データ全点数Nの周波数スペクトルシフト計算を行うことは、そのデータ点数が多いことから処理に時間が掛かる。たとえば、FFT点数1024、デシメーションレート=24とするとNは少なくとも24576となる。また、連続波ドップラモード時には、速度計測処理が200ms毎に実行されるのに対し、この処理に約1秒間のサンプリングデータ列を用いるため、図11(A)に示すように、この1秒間のサンプリングデータ列のうち800ms分は過去の速度計測処理において用いたものとオーバーラップしている。そこで、このように過去において既に周波数スペクトルシフトされたオーバーラップデータについてはそのままそのデータを用い、新たに収集されたデータに対してのみ、上記周波数シフト処理を実行して処理時間を短縮する。ここで、前回の周波数スペクトルシフト処理に用いた離散複素指数関数c(n)の最後の引数nと今回の処理に用いる離散複素指数関数c(n)の先頭の引数nとが連続するようにすれば周波数スペクトルシフトされたデータ列X(n)の位相スペクトルも連続するが、nを無制限に大きくすることができないため、サーキュラバッファ構造のメモリに記憶する。このときFFTに用いるデータ列間で位相のシフトが発生するが、このドップラソナーにおけるFFT処理は周波数のパワースペクトルのピークに基づいてドップラシフト周波数を求めるためのものであり、nが途中でリセットされても〔数3〕から明らかなようにパワースペクトルは保存されているため以後の処理に支障はない。
また、スペクトル帯域を限定するためのフィルタ(s17の処理)としてFIRフィルタが一般的に用いられるが、サンプリングデータx(n)の全てにFIRフィルタリング処理をしても、そのうちデシメーションレートDのデシメーションによって選択される1/Dのデータしか使用されないため無駄である(なお、この説明におけるx(n)は、上記周波数シフトされたX(n)を表すものとする。)。このため、各データ毎にフィルタリング処理をせず、デシメーションによって破棄されるデータをスキップし、D番目のデータ毎にフィルタ処理をするようにする。この破棄されるデータをスキップするFIRフィルタ処理をブロック図で示すと図12のようになる。この図において左端の入力端から1データずつ連続して入力されるサンプリングデータx(n)のうち、デシメーションによって選択されたデータx(Dn)が入力されたときのみ、FIRフィルタが機能し、各タップの演算部はフィルタ係数を乗算したデータを出力する。これらのデータを加算すればフィルタリングされたx(Dn)のデータを得ることができる。デシメーションによって破棄されるデータが入力端から入力された場合には上段の遅延回路(シフトレジスタ)のみ機能し、各タップの演算部の動作を休止させていればよい。これによって、フィルタリングとデシメーションを同時に行うことができる。
ここで、FIRフィルタは、x(n)の実数部および虚数部に対応して2系統が必要であるが、x(n)の各項はそれぞれ実数部または虚数部のみ有意な値をもち他方は0である(図10参照)。したがって、Dが2の倍数であればFIRフィルタから値を出力するときに、有意な実数部と有意な虚数部が与えられるフィルタ係数、および、0である実数部と0である虚数部が与えられるフィルタ係数は決まっている。したがって、Dを2の倍数になるように設定すれば、必ず0が与えられるフィルタ係数およびその演算処理を省略することができる。
さらに、上述したように周波数シフトは、入力されたサンプリングデータを実数部または虚数部に転記するとともにその符号を制御する処理であり、4データ毎に、(+1,−j,−1,+j)または(+1、+j、−1、−j)の処理を繰り返し行うものである。そこで、この実数部または虚数部として取り出す処理および正負符号をフィルタ係数として内蔵することにより、周波数シフト、フィルタリング、デシメーションを同時に処理することができる。
図13に周波数シフト、フィルタリング、デシメーションの一括処理であるFFT前処理プロセスのブロック図を示す。この図は、フィルタ長が奇数であり、且つ、(L+1)/2が偶数の場合を示している。この処理の条件としてデシメーションレートDが周波数シフト処理の繰り返しステップである4の倍数であることが要求される。デシメーションによってD個に1個の割合で選択されるデータx(Dn)が入力されたとき、このx(Dn)からフィルタ長だけ逆上ったデータx(Dn−(L−1))までのデータに対して、フィルタ係数h(0),h(1),−h(2),−h(3),……,−h(L−1)を乗算し、x(Dn)・h(0)+x(Dn−2)・(−h(2))+……+x(Dn(L−1))・(−h(L−1))をフィルタ出力の実数部として出力する。且つ、x(Dn−1)・h(1)+x(Dn−3)・(−h(3))+……+x(Dn(L−1))・h(L−1)をフィルタ出力の虚数部として出力する。これにより、周波数シフトおよびフィルタリングを同時に行うことができる。そして、A/DデータがDだけ進み、次の選択データx(D(n+1))が入力されたとき同様の処理を行う。このD毎の処理によりデシメーションが実行される。このように、A/Dデータを1データずつ進め、デシメーションにより選択されるデータが入力されたとき、フィルタリングおよび周波数シフトを一気に行うようにしたことにより、処理が簡略化されるとともに、無駄な演算処理を全く行うことがないため、演算所要時間を大きく短縮することができる。
なお、上述したようにこの例では、フィルタ長が奇数であり、且つ、(L+1)/2が偶数の場合を示しているが、(L+1)/2が奇数の場合でも、また、フィルタ長が偶数でL/2が奇数・偶数の場合でも同様に処理することが可能である。
なお、図3に示した受信部はアナログ部で一旦中間周波にダウンコンバートしたのちA/D変換する構成になっているが、図14に示すように受波器が受信した信号をそのままA/D変換するようにしてもよい。この場合にはサンプルホールド回路を内蔵した広帯域高速A/D変換器を用いることにこれを実現することができる。この場合でも、サンプリング周波数を適当に選択することにより、実数系統のみでA/D変換でき、且つ、DSP内における周波数スペクトルシフト演算を簡略化することができる。
なお、この実施形態はドップラソナーについてしたが、この発明は、ドップラソナーのみならず、信号を周波数スペクトルに展開し、そのピークを求める処理であればどのような分野にも適用することができる。