JPWO2006126303A1 - 速度計測方法及びこれを用いた速度計測装置 - Google Patents
速度計測方法及びこれを用いた速度計測装置 Download PDFInfo
- Publication number
- JPWO2006126303A1 JPWO2006126303A1 JP2007517722A JP2007517722A JPWO2006126303A1 JP WO2006126303 A1 JPWO2006126303 A1 JP WO2006126303A1 JP 2007517722 A JP2007517722 A JP 2007517722A JP 2007517722 A JP2007517722 A JP 2007517722A JP WO2006126303 A1 JPWO2006126303 A1 JP WO2006126303A1
- Authority
- JP
- Japan
- Prior art keywords
- order
- speed
- expansion coefficient
- complex
- legendre
- 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
Links
Images
Classifications
-
- 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/88—Sonar systems specially adapted for specific applications
- G01S15/89—Sonar systems specially adapted for specific applications for mapping or imaging
- G01S15/8906—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
- G01S15/8979—Combined Doppler and pulse-echo imaging systems
- G01S15/8981—Discriminating between fixed and moving objects or between objects moving at different speeds, e.g. wall clutter filter
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/06—Measuring blood flow
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/48—Diagnostic techniques
- A61B8/488—Diagnostic techniques involving Doppler signals
-
- 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/88—Sonar systems specially adapted for specific applications
- G01S15/89—Sonar systems specially adapted for specific applications for mapping or imaging
- G01S15/8906—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
- G01S15/8977—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using special techniques for image reconstruction, e.g. FFT, geometrical transformations, spatial deconvolution, time deconvolution
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Surgery (AREA)
- Heart & Thoracic Surgery (AREA)
- Veterinary Medicine (AREA)
- Public Health (AREA)
- Biophysics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Acoustics & Sound (AREA)
- Animal Behavior & Ethology (AREA)
- General Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Hematology (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
Ｎ個の時系列信号を０次から（Ｎ−１）次までの離散的ルジャンドル関数を基底として展開するステップ（Ｓ４）と、（２ｎ−１）次の展開係数と（２ｎ＋１）次の展開係数との線形結合に虚数単位を乗じて、２ｎ次の展開係数と線形結合した２ｎ次の複素展開係数と、（２ｎ＋１）次の展開係数に虚数単位を乗じた後、２ｎ次の展開係数と（２ｎ＋２）次の展開係数と線形結合した（２ｎ＋１）次の複素展開係数とを算出するステップ（Ｓ５）と、各複素展開係数のうち最大の絶対値を持つ係数の次数ｍを決定する次数決定ステップ（Ｓ４）と、次数ｍに対応した各展開係数間あるいは各複素展開係数間の自乗和の比から運動反射物に関する符号付き速度信号を算出するステップ（Ｓ８）とからなる速度計測方法及びこの方法を実施する速度計測装置。
Description
本発明は、パルス状の反射波を用いた速度計測方法及びこれを用いた速度計測装置に関する。
パルス・ドップラ速度計では、複数のパルス波を運動反射物に送信し、図１のような受信された複数のエコー信号が解析される。すなわち、各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べた時系列信号が解析され、運動反射物に関する速度などの情報が得られる。最も一般的な信号処理方法としては、これらの受信信号を直交検波し、複素時系列信号として解析する。このようなパルス・ドップラ速度計は、生体内部の血流などを検出・描画する超音波診断装置、雨雲などを検出・描画する気象レーダー、飛行物体を検出する航空レーダーなどとして広く用いられている。
もし、反射物が１つしか存在しなければ、前記の時系列信号の位相回転速度 ΔΦ／Δtすなわち符号つき角周波数から、反射物が送受信器に近づく運動速度ｖを、次式のように容易に求めることができる。
ｖ＝λΔΦ／Δt／２ （１）
ここで、λは波長、Δtはパルス送信の時間間隔、ΔΦは位相回転角である。運動速度ｖの符号は、反射物が近づく場合には、運動速度ｖの符号が正となり、遠ざかる場合には負となる。
もし、反射物が１つしか存在しなければ、前記の時系列信号の位相回転速度 ΔΦ／Δtすなわち符号つき角周波数から、反射物が送受信器に近づく運動速度ｖを、次式のように容易に求めることができる。
ｖ＝λΔΦ／Δt／２ （１）
ここで、λは波長、Δtはパルス送信の時間間隔、ΔΦは位相回転角である。運動速度ｖの符号は、反射物が近づく場合には、運動速度ｖの符号が正となり、遠ざかる場合には負となる。
ところで、前記のような実際の受信エコー中では、静止反射物によるエコー信号、いわゆるクラッタ信号の強度が、着目する運動反射物によるエコー信号の強度よりも、一般に、数桁以上大きい。また、静止反射物によるエコー信号は、実際には、時間軸上で完全に静止しているわけではなく、途中の媒質のゆらぎや、医療用超音波診断装置の場合には、血流以外の静止臓器の遅い動きのために、ドリフトしている。このため、静止反射物によるエコー信号の周波数成分は、位相回転速度ゼロの直流成分だけでなく、位相回転速度がゼロでない低周波成分を有している。
このため、前記時系列信号に単純に (１)式の処理を適用しただけでは、運動反射物の速度を検出することはできない。そこで、実際のパルス・ドップラ速度計は、静止反射物によるエコー信号を抑圧し、運動反射物によるエコー信号を相対的に強調するＭＴＩ(Moving Target Indicator) 処理とよばれる信号処理を行った後に、速度の検出あるいは分析を行うように構成されている。
ＭＴＩ処理として最もよく知られているのは、時間領域において畳み込み積和で表現される通常の低域遮断フィルタを用いる処理である。この処理には次の二つの欠点がある。
１） Ｎ１個の時系列データ点を入力とする低域遮断フィルタを用いる場合、後段の速度検出・分析処理部に入力されるデータ点が (Ｎ１−１) 個目減りする。
２） 急峻な遮断特性のフィルタを得にくい。
ＭＴＩ処理として最もよく知られているのは、時間領域において畳み込み積和で表現される通常の低域遮断フィルタを用いる処理である。この処理には次の二つの欠点がある。
１） Ｎ１個の時系列データ点を入力とする低域遮断フィルタを用いる場合、後段の速度検出・分析処理部に入力されるデータ点が (Ｎ１−１) 個目減りする。
２） 急峻な遮断特性のフィルタを得にくい。
パルス・ドップラ速度計の場合は、運動反射物によるエコー信号を保ちながら、静止反射物のドリフト成分を除去する急峻な遮断特性が必要であり、特に、２）が問題となる。また、Ｎ回の送受信により得られるＮ点の時系列信号から速度検出・分析処理をするとき、データ点が目減りする１）の問題は、実際に速度検出・分析演算に用い得るデータ点が(Ｎ−Ｎ１＋１) 点に減ってしまうことを意味する。これは、医療診断用超音波血流描画装置のように実時間性が重要な応用では、望ましいことではない。
この問題を解決する処理としては、非特許文献１に記載されているようにPolynomial Regression Filter （多項式フィット・フィルタ）が提案されている。これは、時系列信号に、０次式（定数）、１次式、２次式、・・・・、Ｍ次式を順次最小２乗フィットし、フィットした成分を引き去ることにより、元の時系列信号のもつドリフト成分を除去する処理である。この処理は、Ｎ点の入力時系列信号にＮ×Ｎの正方行列を掛けることで表現される。このため、出力信号として同じＮ点の時系列信号が得られ、
１） データ点の目減りがない。
また、その低域遮断特性は、フィットする最大次数Ｍの大きさによるが、
２） 同等の遮断周波数をもつ前記の低域遮断フィルタと比較すると、はるかに急峻な遮断特性をもつ。
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Ｖol. 42., 927-937頁 (1995年)
１） データ点の目減りがない。
また、その低域遮断特性は、フィットする最大次数Ｍの大きさによるが、
２） 同等の遮断周波数をもつ前記の低域遮断フィルタと比較すると、はるかに急峻な遮断特性をもつ。
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Ｖol. 42., 927-937頁 (1995年)
ところが、この多項式フィット・フィルタによる処理の後段において、ｖ＝λΔΦ／Δｔ／２式を用いて運動物体の速度の検出・分析を行うと、遮断周波数の近くで速度演算誤差が、低域遮断フィルタと比較しても大きくなってしまうという問題点がある。このため、多項式フィット・フィルタが有する急峻な遮断特性を十分に使うことができない。
そこで、本発明は、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減することができる速度計測方法及びこれを用いた速度計測装置を提供することを課題とする。
そこで、本発明は、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減することができる速度計測方法及びこれを用いた速度計測装置を提供することを課題とする。
本発明は、速度計測対象物に対してＮ個のパルス波を所定間隔で送信し、受信エコー信号を送信時刻の順番に並べたＮ個の時系列信号を０次から（Ｎ−１）次までの離散的ルジャンドル関数を基底として展開し、この展開係数を用いて速度を算出する。特に、１次以上（Ｎ−１）次以下の（２ｎ−１）次の展開係数と（２ｎ＋１）次の展開係数との線形結合に虚数単位を乗じて、２ｎ次の展開係数と線形結合した２ｎ次の複素展開係数と、（２ｎ＋１）次の展開係数に虚数単位を乗じた後、２ｎ次の展開係数と（２ｎ＋２）次の展開係数と線形結合した（２ｎ＋１）次の複素展開係数とを算出し、各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数ｍを決定し、次数ｍに対応した前記各展開係数間あるいは各複素展開係数間の自乗和の比から運動反射物に関する符号付き速度信号を算出する。
パルス波が照射された運動反射物からは位相回転速度（角周波数）ΔΦ／Δtを有したエコー信号が得られ、ドリフト振動を伴った静止反射物からは低周波成分を有したクラッタ信号が得られる。このクラッタ信号の強度は、エコー信号の強度よりも数桁大きい。エコー信号を次数が大きいほど周波数原点付近で急速にゼロに収束する特性を有する離散的ルジャンドル関数を基底として展開し、各展開係数の横軸を位相回転速度（角周波数）に対応させることによって、運動反射物のエコー信号と静止対象物の低周波成分とが区別される。つまり、静止反射物と運動反射物との反射信号を区別することができる。
本発明によれば、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減できる。
（フィルタの原理）
各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べたＮ点の離散的時系列信号から速度分析を行う場合を考える。このような、パルス・ドップラ速度分析用として典型的な時系列複素信号は、受信されたエコー信号に９０°位相の異なる２つの搬送周波数信号を乗じて得られる１組の信号を、虚数単位を係数として線型結合して複素信号とし、この信号の、各々のパルス送信時刻を基準とする時相の等しい部分を、送信時刻の順番にならべることにより得られる（図１参照）。
このような離散的時系列信号のもつ周波数あるいは位相の分析における常套手段は、離散的フーリエ展開である。具体的には、ｎ＝ １，２，・・・， Ｎ とならんだ時系列信号を、周波数を表す指数ｋ＝−Ｎ／２，・・，０，・・，Ｎ／２の複素正弦波系
Cs(n，k) = cos [ kπ(２ｎ−Ｎ−１)／Ｎ ] + j sin [ kπ(２ｎ−Ｎ−１)/Ｎ ] （２）
で展開する。j は虚数単位である。
負のｋ は反射物が遠ざかる向きの速度に対応し、正のｋ は近づく向きの速度に対応する。ｋ＝±Ｎ／２ はナイキスト限界であり、正負が区別されないので、(２)式により表される (Ｎ＋１)個の複素正弦波関数のうち独立なものはＮ個となり、これらは、直交関数系を形成している。
各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べたＮ点の離散的時系列信号から速度分析を行う場合を考える。このような、パルス・ドップラ速度分析用として典型的な時系列複素信号は、受信されたエコー信号に９０°位相の異なる２つの搬送周波数信号を乗じて得られる１組の信号を、虚数単位を係数として線型結合して複素信号とし、この信号の、各々のパルス送信時刻を基準とする時相の等しい部分を、送信時刻の順番にならべることにより得られる（図１参照）。
このような離散的時系列信号のもつ周波数あるいは位相の分析における常套手段は、離散的フーリエ展開である。具体的には、ｎ＝ １，２，・・・， Ｎ とならんだ時系列信号を、周波数を表す指数ｋ＝−Ｎ／２，・・，０，・・，Ｎ／２の複素正弦波系
Cs(n，k) = cos [ kπ(２ｎ−Ｎ−１)／Ｎ ] + j sin [ kπ(２ｎ−Ｎ−１)/Ｎ ] （２）
で展開する。j は虚数単位である。
負のｋ は反射物が遠ざかる向きの速度に対応し、正のｋ は近づく向きの速度に対応する。ｋ＝±Ｎ／２ はナイキスト限界であり、正負が区別されないので、(２)式により表される (Ｎ＋１)個の複素正弦波関数のうち独立なものはＮ個となり、これらは、直交関数系を形成している。
反射物が１つしか存在せず、それがナイキスト限界範囲内の一定速度で運動している場合を想定し、運動速度の関数として、前記展開により得られる離散的フーリエ展開係数を計算し、Ｎ＝８の場合について、その絶対値を図２にプロットした。横軸は、運動速度に対応する位相回転速度を、パルス送信の時間間隔の逆数 Pulse Repetition Frequency (PRF) の２π倍を単位として表示した。図２では展開係数を表す８個の関数が重なって見づらいので、そのうち４個ずつを見分けやすい線により図３及び図４に示した。
送信時間間隔のＮ倍の時間で、位相が整数回回転するような運動速度では、ゼロでない展開係数はただ１つであり、他の展開係数はすべてゼロとなり、きれいに速度分析されている。ところが、それ以外の位相回転が整数回とならない一般の運動速度では、それぞれの展開係数が、横軸上ピークからはるかに離れた運動速度においても、ピークの１／１０程度の大きさの絶対値をもっている。これは、速度分析計として−２０ｄＢ程度の遠隔的クロストークをもつことを意味する。検出・分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号がドリフトするとき、これはクラッタ信号と運動反射物によるエコー信号とを区別できないという重大な問題を生じさせる。
この問題は、Ｎ個の時系列信号に、ハニング関数のような緩やかに立ち上がって、下がる窓関数による重みをかけることにより、かなり抑圧できるが、本質的には解決しない。また、この方法は、時系列信号の数が大幅に目減りしたのと同様の、望ましくない効果を生じ、特に、時系列信号の数が元々少ない場合には適さない方法である。
具体例として、ナイキスト限界の２/３の一定速度で運動する検出・分析対象である運動反射物があり、その１０００倍のエコー信号振幅をもつクラッタ信号が、ナイキスト限界速度の１／１００でドリフトしている場合を考える。図５には、このときの前記速度分析計の出力である速度スペクトルを示す。検出・分析対象である運動反射物のスペクトルは、速度ゼロ付近にピークをもつクラッタ信号からのクロストーク成分に覆いかくされてしまい、この出力スペクトルの中では検出不可能である。このことが、まさに、クラッタ信号を遮断するＭＴＩフィルタが必要な最も本質的な理由である。
次に、ＭＴＩ処理として最もよく知られ、時間領域において畳み込み積和で表現される通常の低域遮断フィルタについて、その通過振幅特性の典型例を、図６中に実線で示す。このようなフィルタにｎ１ = １,・・・, （Ｎ＋Ｍ−１） とならんだ時系列信号Ｓ１(ｎ１) を、入力すると、出力としてｎ０＝１,・・・, Ｎ とならんだ時系列信号Ｓ０(ｎ０)が得られる。
Ｓ０(ｎ０)＝ΣＦ（ｍ）Ｓ１(ｎ０＋ｍ−１) (３)
ここで、Σは ｍ＝１，・・・，Ｍの和を表す。図６では、低域遮断フィルタとして最も簡単なＭ＝３，Ｆ(１)＝−１， Ｆ（２）＝２, Ｆ（３）＝−１の場合について実線で示した。図５と同様の入力信号を、この低域遮断フィルタを通した後、図５と同様に速度分析すると、図７中の実線のようになる。低域遮断フィルタの働きによりクラッタ信号の振幅が１／２０００倍程度に抑圧された結果、ナイキスト限界の０．３倍から０．４倍の速度にピークをもつ運動反射物のスペクトルが見えるようになっている。
Ｓ０(ｎ０)＝ΣＦ（ｍ）Ｓ１(ｎ０＋ｍ−１) (３)
ここで、Σは ｍ＝１，・・・，Ｍの和を表す。図６では、低域遮断フィルタとして最も簡単なＭ＝３，Ｆ(１)＝−１， Ｆ（２）＝２, Ｆ（３）＝−１の場合について実線で示した。図５と同様の入力信号を、この低域遮断フィルタを通した後、図５と同様に速度分析すると、図７中の実線のようになる。低域遮断フィルタの働きによりクラッタ信号の振幅が１／２０００倍程度に抑圧された結果、ナイキスト限界の０．３倍から０．４倍の速度にピークをもつ運動反射物のスペクトルが見えるようになっている。
同じ目的のＭＴＩフィルタとして、前記の多項式フィット・フィルタを設計したときの通過振幅特性の典型例を、図６中に破線で示す。この例では、０次から３次式までをフィットして引き去った。前記の畳み込み型フィルタよりも急峻な低域遮断特性が見てとれる。また、図５と同様の入力信号を、この低域遮断フィルタを通した後、図５と同様に速度分析した結果を、図７中に破線で示した。低域遮断特性の急峻さのために、クラッタ信号が完全に抑圧されている。多項式フィット・フィルタ特有の急峻な低域遮断特性は以下のように理解することができる。
連続座標系における狭義ルジャンドル多項式については、岩波数学公式III特殊関数８２−８５頁 に記述されており、例えば、８次まで記すと
Ｐ0（ｘ）＝１
Ｐ1（ｘ）＝ｘ＝cosθ
Ｐ2（ｘ）＝（１／２）（３ｘ2−１）＝（１／４）（３cos２θ＋１）
Ｐ3（ｘ）＝（１／２）（５ｘ3−３ｘ）＝（１／８）（５cos３θ＋３cosθ）
Ｐ4（ｘ）＝（１／８）（３５ｘ4−３０ｘ2＋３）＝（１／６４）（３５cos４θ＋２０cos２θ＋９）
Ｐ5（ｘ）＝（１／８）（６３ｘ5−７０ｘ3＋１５ｘ）＝（１／１２８）（６３cos５θ＋３５cos３θ＋３０cosθ）
Ｐ6（ｘ）＝（１／１６）（２３１ｘ6−３１５ｘ4＋１０５ｘ2−５）＝（１／５１２）（２３１cos６θ＋１２６cos４θ＋１０５cos２θ＋５０）
Ｐ7（ｘ）＝（１／１６）（４２９ｘ7−６３９ｘ5＋３１５ｘ3−３５ｘ）＝（１／２10）（４２９cos７θ＋２３１cos５θ＋１８９cos３θ＋１７５cosθ）
Ｐ8（ｘ）＝（１／１２８）（６４３５ｘ8−１２０１２ｘ6＋６９３０ｘ4−１２６０ｘ2＋３５）＝（１／２14）（６４３５cos８θ＋３４３２cos６θ＋２７７２cos４θ＋２５２０cos２θ＋１２２５）
である。離散的座標系についても使えるように、一般化すると、偶数（２ｎ）次ルジャンドル多項式は、２ｎより低い次数のルジャンドル多項式すべてと直交する最高次数２ｎの偶関数と定義することができ、奇数（２ｎ＋１）次ルジャンドル多項式は、(２ｎ＋１)より低い次数のルジャンドル多項式すべてと直交する最高次数(２ｎ＋１)の奇関数と定義することができる。
Ｐ0（ｘ）＝１
Ｐ1（ｘ）＝ｘ＝cosθ
Ｐ2（ｘ）＝（１／２）（３ｘ2−１）＝（１／４）（３cos２θ＋１）
Ｐ3（ｘ）＝（１／２）（５ｘ3−３ｘ）＝（１／８）（５cos３θ＋３cosθ）
Ｐ4（ｘ）＝（１／８）（３５ｘ4−３０ｘ2＋３）＝（１／６４）（３５cos４θ＋２０cos２θ＋９）
Ｐ5（ｘ）＝（１／８）（６３ｘ5−７０ｘ3＋１５ｘ）＝（１／１２８）（６３cos５θ＋３５cos３θ＋３０cosθ）
Ｐ6（ｘ）＝（１／１６）（２３１ｘ6−３１５ｘ4＋１０５ｘ2−５）＝（１／５１２）（２３１cos６θ＋１２６cos４θ＋１０５cos２θ＋５０）
Ｐ7（ｘ）＝（１／１６）（４２９ｘ7−６３９ｘ5＋３１５ｘ3−３５ｘ）＝（１／２10）（４２９cos７θ＋２３１cos５θ＋１８９cos３θ＋１７５cosθ）
Ｐ8（ｘ）＝（１／１２８）（６４３５ｘ8−１２０１２ｘ6＋６９３０ｘ4−１２６０ｘ2＋３５）＝（１／２14）（６４３５cos８θ＋３４３２cos６θ＋２７７２cos４θ＋２５２０cos２θ＋１２２５）
である。離散的座標系についても使えるように、一般化すると、偶数（２ｎ）次ルジャンドル多項式は、２ｎより低い次数のルジャンドル多項式すべてと直交する最高次数２ｎの偶関数と定義することができ、奇数（２ｎ＋１）次ルジャンドル多項式は、(２ｎ＋１)より低い次数のルジャンドル多項式すべてと直交する最高次数(２ｎ＋１)の奇関数と定義することができる。
例として、一連の８点からなる時系列信号について、０次から７次までの離散的ルジャンドル関数を計算し、偶数次を図８に、奇数次を図９に示した。０次・１次を実線、２次・３次を破線、４次・５次を１点鎖線、６次・７次を点線にて示した。振幅は、ＲＭＳ（２乗平均の平方根）値が１となるよう規格化した。符号は、最低次項の符号が正となるよう、すなわち、偶数次多項式については定数項、奇数次多項式については１次項の符号が正となるよう定めた。図１０には、その周波数スペクトルをプロットした。周波数成分の振幅の絶対値を、個々のルジャンドル展開係数の最大値で規格化して示した。図１０では８個の関数が重なってしまって見づらいので、そのうち４個ずつを見分けやすい線により図１１及び図１２に示した。Ａ（０），Ａ（７）を実線、Ａ（１），Ａ（６）を破線、Ａ（２），Ａ（５）を1点鎖線、Ａ（３），Ａ（４）を点線にて示した。
図１０から図１２を、図２から図４と同様、入力複素正弦波に対する周波数応答として眺めてみる。ピークを主応答、それ以外を副応答とみる。ルジャンドル関数は、図２から図４のフーリエ展開系に比べ、主応答より高周波側の副応答は大きいが、低周波側には副応答をもたず、次数が大きいほど周波数原点付近で急速にゼロに収束する。ここでは解析的証明は省略するが、周波数の関数としてのｍ次ルジャンドル関数は、周波数原点において０次から(ｍ−１)次導関数まですべてが０となる。
前記の多項式フィット・フィルタは、入力信号を、この離散的ルジャンドル関数を基底として展開し、０からＭ次までを除き、(Ｍ＋１)次以上の成分の和を出力信号とするフィルタであるということができる。したがって、出力信号は、周波数原点付近で周波数の(Ｍ＋１)乗に比例する挙動でゼロに収束していく。図６及び図７中に破線で示された例のように、多項式フィット・フィルタの低域遮断特性が優れているのは、このようなルジャンドル関数の本質的な特性に由来するものである。このように考えると、クラッタという周波数原点付近に桁違いに大きな振幅をもつ妨害成分の影響を抑えて、着目する運動物体を検出したり、その速度を分析したりしなければならないドップラ速度検出には、多項式フィット・フィルタやルジャンドル関数展開を利用する方式は、本質的に適しているといえる。
フーリエ展開系では、複素正弦波を基底とすることにより、図２から図５や図７に示したように符号つきの速度分析が可能であることは広く知られており、そのような分析方法が広く用いられている。ところが、ルジャンドル展開系は、前記のように、ドップラ速度検出に極めて適した性質をもつものの、このまま用いたのでは、図１０から図２０より見てとれるように、絶対値が等しいが符号の異なる速度に対して等しい応答を示すので、符号つきの速度分析は不可能である。
そこで、次数が１つだけ異なる偶数次と奇数次のルジャンドル関数を虚数単位 j を係数として線型結合して、複素ルジャンドル関数を形成し、これにより符号つきの速度分析を可能とする。これは、偶数次のルジャンドル関数が余弦波に、奇数次のルジャンドル関数が正弦波に似た変化を示すことに着目し、余弦関数と正弦関数を虚数単位 j を係数として線型結合することにより、位相角をｊ倍した指数関数、すなわち複素正弦関数が得られ、この位相角の増減から符号つき速度が得られるのと同様に、複素ルジャンドル関数から符号つき速度を得ようとするものである。
具体的には、一連のＮ点からなる時系列信号を、まず、０次から (Ｎ−１)次の離散的ルジャンドル関数を基底として展開して、Ａ（０）からＡ（Ｎ−１）までの展開係数を得て、これをもとに、一連の複素ルジャンドル係数
Ｃ（±(２ｎ−３／２)）＝Ａ（２ｎ−２）±ｊＡ（２ｎ−１） (１≦２ｎ−１≦Ｎ-1)
Ｃ（±(２ｎ−１／２)）＝Ａ（２ｎ）±j Ａ（２ｎ−１） (２≦２ｎ≦Ｎ−1;ｎ:自然数, 複号同順) (４)
を算出し、この係数の相対的な大きさにより符号つき速度分析を行う。
具体的には、一連のＮ点からなる時系列信号を、まず、０次から (Ｎ−１)次の離散的ルジャンドル関数を基底として展開して、Ａ（０）からＡ（Ｎ−１）までの展開係数を得て、これをもとに、一連の複素ルジャンドル係数
Ｃ（±(２ｎ−３／２)）＝Ａ（２ｎ−２）±ｊＡ（２ｎ−１） (１≦２ｎ−１≦Ｎ-1)
Ｃ（±(２ｎ−１／２)）＝Ａ（２ｎ）±j Ａ（２ｎ−１） (２≦２ｎ≦Ｎ−1;ｎ:自然数, 複号同順) (４)
を算出し、この係数の相対的な大きさにより符号つき速度分析を行う。
図１３から図１５には、図１０から図１２と同様Ｎ＝８の場合について、複素正弦波入力に対する出力としての複素ルジャンドル係数の絶対値を、符号つき速度に対応する位相回転速度の関数として示した。複素ルジャンドル係数の絶対値は、それぞれの最大値で規格化してプロットした。図１３中の実線はＣ（６．５）、破線はＣ（−６．５）、図１４中の実線はＣ（５．５）、破線はＣ（−５．５）、図１５中の実線はＣ（４．５）、破線はＣ（−４．５）である。これらの図から、複素ルジャンドル係数が、速度に対し正負を区別して応答していることがわかる。なお、このような速度に対する応答をもつ複素ルジャンドル系がＮ≦３５の範囲で構築できることが確認されている。
図１３から図１５に示した複素ルジャンドル係数の応答は、ルジャンドル展開係数の特徴である位相回転速度原点すなわちドップラ周波数原点のまわりの極めて優れたカットオフ特性を示している反面、図２から図４に示したフーリエ展開係数の応答と比較して、周波数サイドローブレベルが高めである。すなわち、周波数サイドローブレベルは、フーリエ展開係数の図２から図４の例では、０．２程度であるのに対し、複素ルジャンドル係数の図１４の例では、０．６近くに達している。周波数サイドローブレベルが高いことは、S/Ｎ比すなわちノイズレベルに対する信号レベルが低い場合に、速度検出誤差を生む可能性が高いことを意味するので、周波数サイドローブレベルは抑圧することが望ましい。
そこで、Ａ（２ｎ）をもとにＣ（±２ｎ）を生成する際、Ａ（２ｎ）と同じ周波数において最大値をとるようＡ（２ｎ−１）とＡ（２ｎ−１）とを線型結合して用い、また、Ａ（２ｎ−１）をもとにＣ（±(２ｎ−１)）を生成する際、Ａ（２ｎ−１）と同じ周波数において最大値をとるようＡ（２ｎ−２）とＡ（２ｎ）とを線型結合して用いる。
すなわち、
Ｃ（±(２ｎ−１)）＝α(２ｎ−１)×Ａ（２ｎ−２) ＋β(２ｎ−１)×Ａ（２ｎ）±j Ａ（２ｎ−１） (１≦２ｎ−１≦Ｎ−１)
Ｃ（±２ｎ）＝Ａ（２ｎ）±j [α(２ｎ)×Ａ（２ｎ−１）+ β(２ｎ)×Ａ（２ｎ＋１） ] (２≦２ｎ≦Ｎ−１；ｎ：自然数，複号同順；α及びβ：正の実数 ) (５)
となる。
すなわち、
Ｃ（±(２ｎ−１)）＝α(２ｎ−１)×Ａ（２ｎ−２) ＋β(２ｎ−１)×Ａ（２ｎ）±j Ａ（２ｎ−１） (１≦２ｎ−１≦Ｎ−１)
Ｃ（±２ｎ）＝Ａ（２ｎ）±j [α(２ｎ)×Ａ（２ｎ−１）+ β(２ｎ)×Ａ（２ｎ＋１） ] (２≦２ｎ≦Ｎ−１；ｎ：自然数，複号同順；α及びβ：正の実数 ) (５)
となる。
図１６から図１７には、図１３から図１５と同様Ｎ＝８の場合について、ルジャンドル係数Ａ（４）、Ａ（５）、Ａ（６）から複素ルジャンドル係数Ｃ（±５）を生成する例を示した。図１６中の実線はＡ（５）の絶対値のドップラ周波数応答である。一方、点線は、そのＡ（５）と同じ周波数において最大値をとるようにＡ（４）とＡ（６）とを線型結合して生成した係数の絶対値のドップラ周波数応答である。図１７は、そのＡ（４）とＡ（６）との線型結合を実部、Ａ（５）を虚部として得られる複素ルジャンドル係数Ｃ（±５）の絶対値のドップラ周波数応答である。図中実線はＣ（５）であり、破線はＣ（−５）である。周波数サイドローブレベルは、０．４程度に抑えられている。これは、図１４に示されたＣ（±５．５）や図１５に示されたＣ（±４．５）の周波数サイドローブレベル０．６程度に比べ、大きく改善されている。同様に生成した複素ルジャンドル係数Ｃ（±６）、Ｃ（±４）、Ｃ（±３）の絶対値のドップラ周波数応答を図１８から図２０に示した。これらは図１７と同様に、実線がＣ（６）、Ｃ（４）、Ｃ（３）であり、破線がＣ（−６）、Ｃ（−４）、Ｃ（−３）である。周波数サイドローブレベルは、やはり０．４程度に抑えられている。すなわち、Ａ（２ｎ−１）と同じ角周波数において最大値をとるようＡ（２ｎ−２）とＡ（２ｎ）とを線型結合することによって、サイドローブレベルが０．６から０．４程度に抑えられる。
図２１には、ｍを自然数とするとき、(ｍ−１)次から(Ｎ−１)次までのルジャンドル係数の自乗和とｍ次から(Ｎ−１)次までのルジャンドル係数の自乗和との比
Ｒ（ｍ）＝[Ａ（ｍ）2＋・・・＋Ａ（Ｎ−１)2] ／ [Ａ（ｍ−１）2＋Ａ（ｍ）2＋・・・＋Ａ（Ｎ−１)2] (１≦ｍ≦Ｎ−１) (６)
を、前記の例と同様Ｎ＝８の場合について示した。これらの曲線は、ドップラ周波数原点から始まり、原点から離れるにともない単調増加して、一旦最大値をとり、その後リップルする。その単調増加部分を用いることにより、図２２に示すような、速度校正曲線を得ることができる。絶対値が最大の複素ルジャンドル係数の次数が求まると、その次数に対応する速度校正曲線を用いて、ルジャンドル係数の自乗和の比からドップラ速度を算出すればよい。
Ｒ（ｍ）＝[Ａ（ｍ）2＋・・・＋Ａ（Ｎ−１)2] ／ [Ａ（ｍ−１）2＋Ａ（ｍ）2＋・・・＋Ａ（Ｎ−１)2] (１≦ｍ≦Ｎ−１) (６)
を、前記の例と同様Ｎ＝８の場合について示した。これらの曲線は、ドップラ周波数原点から始まり、原点から離れるにともない単調増加して、一旦最大値をとり、その後リップルする。その単調増加部分を用いることにより、図２２に示すような、速度校正曲線を得ることができる。絶対値が最大の複素ルジャンドル係数の次数が求まると、その次数に対応する速度校正曲線を用いて、ルジャンドル係数の自乗和の比からドップラ速度を算出すればよい。
検出・分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号が存在し、それがドリフトするときには、そのドリフトの程度が大きくなるに応じて、Ｃ（±１）、Ｃ（±２）、……と、次数の小さな係数を順次除いた複素ルジャンドル係数系を用いて、速度検出・分析を行えばよい。このようにルジャンドル係数の遮断次数を制御することにより、クラッタ信号の影響を抑圧しながら、運動反射物の速度検出・分析を行うことができる。
（第１実施形態）
以下、本発明の実施形態について説明する。
図２３は、本発明の一実施形態である速度計測装置を使用した超音波診断装置のブロック図であり、血流描画機能を備えている。超音波探触子１を構成する各素子は、切り替えスイッチ群２を介して、送波ビームフォーマ３と受波ビームフォーマ１０に接続されている。送波ビームフォーマ３では、送受信シークエンス制御部６による制御に従って、送信波形メモリ５から送信波形選択部４により選択されて読み出された波形を用い、各素子を通じて送信されたときに指向性を持つ超音波パルスとなるような信号が生成される。この信号が、超音波探触子１の各素子により超音波パルスに変換されて生体に送信される。生体中で反射あるいは散乱されて超音波探触子１に戻ってきた超音波エコー信号は、各素子に受信されて、電気信号に変換される。受波ビームフォーマ１０では、送受信シークエンス制御部６による制御に従って、指向性を持つ受信感度を生成すべく、各受波信号に遅延時間を与えて互いに加算する。遅延加算により得られた時系列信号は、やはり送受信シークエンス制御部６による制御に従って、受波メモリ選択部１１により選択された受波信号メモリ１２中の１つのバンクへ一旦書き込まれ、ドップラ信号分析すべきＮ個の時系列信号がそろったのちに読み出されて、位相回転検出器１３，１５、ダウン・ミキシング器１４、血流信号検出・分析部１６において速度を検出・分析すべく信号処理される。
以下、本発明の実施形態について説明する。
図２３は、本発明の一実施形態である速度計測装置を使用した超音波診断装置のブロック図であり、血流描画機能を備えている。超音波探触子１を構成する各素子は、切り替えスイッチ群２を介して、送波ビームフォーマ３と受波ビームフォーマ１０に接続されている。送波ビームフォーマ３では、送受信シークエンス制御部６による制御に従って、送信波形メモリ５から送信波形選択部４により選択されて読み出された波形を用い、各素子を通じて送信されたときに指向性を持つ超音波パルスとなるような信号が生成される。この信号が、超音波探触子１の各素子により超音波パルスに変換されて生体に送信される。生体中で反射あるいは散乱されて超音波探触子１に戻ってきた超音波エコー信号は、各素子に受信されて、電気信号に変換される。受波ビームフォーマ１０では、送受信シークエンス制御部６による制御に従って、指向性を持つ受信感度を生成すべく、各受波信号に遅延時間を与えて互いに加算する。遅延加算により得られた時系列信号は、やはり送受信シークエンス制御部６による制御に従って、受波メモリ選択部１１により選択された受波信号メモリ１２中の１つのバンクへ一旦書き込まれ、ドップラ信号分析すべきＮ個の時系列信号がそろったのちに読み出されて、位相回転検出器１３，１５、ダウン・ミキシング器１４、血流信号検出・分析部１６において速度を検出・分析すべく信号処理される。
読み出されたＮ個の時系列受信信号は、各々を得るために行われたパルス送信の時刻を基準とする時相の等しい部分を、送信時刻の順番に並べたのちに処理される。この速度検出アルゴリズムについて、図２４のフローチャート及び図２３を参照して説明する。まず、位相回転検出器１３では、クラッタ信号を含んだままの時系列信号から位相回転を求める。位相回転を検出する最も典型的な方法はｎ＝１,・・・, Ｎとならんだ隣接する時系列複素信号 Ｓ（ｎ），Ｓ（ｎ＋１）をもとに、
Ｐ（ｎ）＝Ｓ（ｎ＋１）Ｓ（ｎ）* ／||Ｓ（ｎ＋１）||／||Ｓ（ｎ）|| (ｎ＝１,・・・,Ｎ−１) (７)
を計算し、Ｐ（ｎ）の平均値Ｐａの位相から平均位相回転速度を求める。
ここで、Ｓ（ｎ）*はＳ（ｎ）の複素共役、||Ｓ（ｎ）|| はＳ（ｎ）の絶対値である。一般に、クラッタ・エコー信号は血流エコー信号より桁違いに信号振幅が大きいので、Ｐａは、概ねクラッタ・エコー信号の位相回転平均値と考えてよい。すなわち、隣接サンプル間位相差平均値よりクラッタ平均速度が計算される（図２４のＳ１）。また、時系列複素信号Ｓ（ｎ）は、受信されたエコー信号に９０°位相の異なる２つの搬送周波数信号を乗じて得られる。
Ｐ（ｎ）＝Ｓ（ｎ＋１）Ｓ（ｎ）* ／||Ｓ（ｎ＋１）||／||Ｓ（ｎ）|| (ｎ＝１,・・・,Ｎ−１) (７)
を計算し、Ｐ（ｎ）の平均値Ｐａの位相から平均位相回転速度を求める。
ここで、Ｓ（ｎ）*はＳ（ｎ）の複素共役、||Ｓ（ｎ）|| はＳ（ｎ）の絶対値である。一般に、クラッタ・エコー信号は血流エコー信号より桁違いに信号振幅が大きいので、Ｐａは、概ねクラッタ・エコー信号の位相回転平均値と考えてよい。すなわち、隣接サンプル間位相差平均値よりクラッタ平均速度が計算される（図２４のＳ１）。また、時系列複素信号Ｓ（ｎ）は、受信されたエコー信号に９０°位相の異なる２つの搬送周波数信号を乗じて得られる。
ダウン・ミキシング器１４では、ここで求めた平均値Ｐａ又はＰａの周辺領域における空間平均値を用い、エコー信号の位相回転平均値がゼロとなるようにミキシング処理が行われる。いいかえれば、ダウン・ミキシングによりクラッタ平均速度による位相回転が打ち消される（図２４のＳ２参照）。すなわち、Ｓ（ｎ）とＰａより
Ｓｄ（ｎ）＝Ｓ（ｎ）Ｐａ* n (ｎ＝１,・・・, Ｎ) (８)
を得る。このダウン・ミキシング処理を行うことにより、クラッタ信号を抑圧する後段の処理をより効果的に行うことができる。
Ｓｄ（ｎ）＝Ｓ（ｎ）Ｐａ* n (ｎ＝１,・・・, Ｎ) (８)
を得る。このダウン・ミキシング処理を行うことにより、クラッタ信号を抑圧する後段の処理をより効果的に行うことができる。
得られた隣接サンプルである時系列複素信号Ｓｄ（ｎ），Ｓｄ（ｎ＋１）をもとに、位相回転検出器１５では、位相回転検出器１３と同様に、
Ｐｄ（ｎ）＝Ｓｄ（ｎ＋１）Ｓｄ（ｎ）* ／ ||Ｓｄ（ｎ＋１）|| ／||Ｓｄ（ｎ）|| (ｎ＝１,・・・, Ｎ−１) (９)
を計算する。この隣接サンプル間位相差最大値より、消し残したクラッタ速度の最大値が評価される（図２４のＳ３）。この位相回転最大値の大きさに応じて、血流信号検出・分析部１６の遮断特性を制御するために、ルジャンドル展開係数の最低次数が決定される（図２４のＳ４）。すなわち、位相回転最大値が大きいときには、小さいときに比べ、ルジャンドル展開係数の遮断次数Ｍがより高く設定され、クラッタ成分が効果的に抑圧される。
Ｐｄ（ｎ）＝Ｓｄ（ｎ＋１）Ｓｄ（ｎ）* ／ ||Ｓｄ（ｎ＋１）|| ／||Ｓｄ（ｎ）|| (ｎ＝１,・・・, Ｎ−１) (９)
を計算する。この隣接サンプル間位相差最大値より、消し残したクラッタ速度の最大値が評価される（図２４のＳ３）。この位相回転最大値の大きさに応じて、血流信号検出・分析部１６の遮断特性を制御するために、ルジャンドル展開係数の最低次数が決定される（図２４のＳ４）。すなわち、位相回転最大値が大きいときには、小さいときに比べ、ルジャンドル展開係数の遮断次数Ｍがより高く設定され、クラッタ成分が効果的に抑圧される。
ここで、本実施形態の特徴部である血流信号検出・分析部１６の動作について、もう少し詳しく述べる。入力された時系列信号 Ｓｄ(１) ,・・・, Ｓｄ（Ｎ） は、まず、図８及び図９のような０から(Ｎ−１)次の離散的ルジャンドル関数に展開され、その展開係数Ａ(０) ,・・・, Ａ(Ｎ−１) が求められる。この計算は、次のような行列演算により容易に行うことができる。ｎ次ルジャンドル関数を行ベクトル Ｌ(ｎ) で表し、これをｎ＝０から(Ｎ−１)までならべてＮ×Ｎ行列ＬＬを生成すると、Ｓｄ(１) ,・・・, Ｓｄ(Ｎ)を要素とする列ベクトルＳｄからＡ(０) ,・・・,Ａ(Ｎ−１)を要素とする列ベクトルＡを得るための行列は、
Ｆ＝(Ｌ・tＬＬ)-1ＬＬ （１０）
と求めることができる。ここで、tＬＬはＬＬの転置行列、ＬＬ-1はＬＬの逆行列を表す。この行列を予め準備しておけば、行列演算
Ａ＝Ｆ・Ｓｄ (１１)
によって迅速に、ルジャンドル展開係数を求めることができる。この展開係数Ａ（０）,・・・, Ａ（Ｎ−１）から、式(５)に従って、Ｃ（±１) ,・・・, Ｃ（±(Ｎ−１))の複素ルジャンドル係数が生成される（図２４のＳ５）。いいかえれば、ｎ次ルジャンドル関数を行ベクトルＬ（ｎ）で表し、これをｎ＝０から（Ｎ−１）まで列方向に並べたＮ×Ｎ行列ＬＬを転置した転置行列ｔＬＬを生成する過程と、転置行列ｔＬＬにＮ×Ｎ行列ＬＬを前から乗算したものの逆行列（ＬＬ・ｔＬＬ)-1を算出する過程と、Ｎ×Ｎ行列ＬＬに逆行列（ＬＬ・ｔＬＬ)-1を前から乗算する過程とによって算出された行列ＦＦを要素が時系列信号である列ベクトルＳｄ（ｎ）に前から乗算して、ルジャンドル展開係数を求めることができる。前記したように別途、位相回転最大値に応じてルジャンドル展開係数の最低次数（遮断次数Ｍ）が決定され（図２４のＳ４）、Ｃ（±１) から Ｃ（±(Ｍ＋１)) までの低次数複素ルジャンドル係数が棄却される。
Ｆ＝(Ｌ・tＬＬ)-1ＬＬ （１０）
と求めることができる。ここで、tＬＬはＬＬの転置行列、ＬＬ-1はＬＬの逆行列を表す。この行列を予め準備しておけば、行列演算
Ａ＝Ｆ・Ｓｄ (１１)
によって迅速に、ルジャンドル展開係数を求めることができる。この展開係数Ａ（０）,・・・, Ａ（Ｎ−１）から、式(５)に従って、Ｃ（±１) ,・・・, Ｃ（±(Ｎ−１))の複素ルジャンドル係数が生成される（図２４のＳ５）。いいかえれば、ｎ次ルジャンドル関数を行ベクトルＬ（ｎ）で表し、これをｎ＝０から（Ｎ−１）まで列方向に並べたＮ×Ｎ行列ＬＬを転置した転置行列ｔＬＬを生成する過程と、転置行列ｔＬＬにＮ×Ｎ行列ＬＬを前から乗算したものの逆行列（ＬＬ・ｔＬＬ)-1を算出する過程と、Ｎ×Ｎ行列ＬＬに逆行列（ＬＬ・ｔＬＬ)-1を前から乗算する過程とによって算出された行列ＦＦを要素が時系列信号である列ベクトルＳｄ（ｎ）に前から乗算して、ルジャンドル展開係数を求めることができる。前記したように別途、位相回転最大値に応じてルジャンドル展開係数の最低次数（遮断次数Ｍ）が決定され（図２４のＳ４）、Ｃ（±１) から Ｃ（±(Ｍ＋１)) までの低次数複素ルジャンドル係数が棄却される。
次に、棄却されないで残された正次数の複素ルジャンドル係数の自乗和と負次数の複素ルジャンドル係数の自乗和とを比較して、速度符号が決定される（図２４のＳ６参照）。具体的に、大きな方の次数の符号を速度の符号とする。さらに、速度符号次数の複素ルジャンドル係数のうち、絶対値最大の係数が検出される（図２４のＳ７参照）。具体的に、正なら正の次数、負なら負の次数の複素ルジャンドル係数と、最大次数(Ｎ−１)のルジャンドル係数のうち、絶対値最大の係数の次数ｍが検出される。
そして、検出したルジャンドル係数に対応する校正曲線を用いて血流速度を算出する（図２４のＳ８参照）。すなわち、次数ｍに対応する図２２に示したような速度校正曲線を用いて、式(６)のように(ｍ−１)次から(Ｎ−１)次までのルジャンドル係数の自乗和とｍ次から(Ｎ−１)次までのルジャンドル係数の自乗和との比からドップラ速度を求める。
得られた血流速度信号は、エコー振幅検出器１７及びエコー振幅圧縮器１８により得られる静止臓器からのエコー信号と共に、スキャンコンバータ１９に入力される。スキャンコンバータでは、入力された複数の信号を適宜重畳して表示器２０にて２次元ないし３次元表示すべく、信号の生成・制御を行う。
次に、本実施形態のドップラ速度計測装置の動作例を従来方式と比較して示す。クラッタ・エコー信号の振幅は、血流エコーの３００倍とし、クラッタ速度については、初期値がゼロで一定の加速度で立ち上がる場合を想定した。時系列信号数Ｎ＝８の場合について、数値計算シミュレーションを行って比較した。従来方式としては、ＭＴＩフィルタとして図６の場合の低域遮断フィルタを用い、その後段において、式(５)を用いて位相回転量平均値Ｐａを得たのと同様の方法で、速度を算出する処理を血流信号検出・分析部１６において行う例を選んだ。さらに、ＭＴＩフィルタを、やはり図５に例を示したような多項式フィット・フィルタに置き換えた場合とも、比較した。この場合、遮断次数Mは、本発明の方法と同様、位相回転検出器１５の出力信号に応じて制御した。
図２５から図２７には、クラッタの最終到達速度がナイキスト限界速度の０．８％である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、及び多項式フィット・フィルタ方式による速度検出結果を示した。横軸に入力としての血流速度、縦軸に出力として検出された速度をとり、有効速度検出範囲を実線、範囲外を点線によりプロットした。また、理想的場合を一点鎖線により示した。このように、クラッタ速度が低い場合には、いずれの方式も、ほぼ正確な速度分析結果を与えているが、本実施形態の方式による誤差が最も小さい。
図２８から図３０には、同様に、クラッタの最終到達速度がナイキスト限界速度の３％である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。このクラッタ速度域では、低域遮断フィルタ方式による速度検出は、ほぼ破綻に至っている。また、多項式フィット・フィルタ方式は、有効速度検出範囲における誤差は大きくないものの、血流速度がそれより小さい場合の誤差が大きく、何らかの対策を加えない限り、使うことができない。これに対し、本実施形態の方式は、有効速度検出範囲において正確な速度検出結果を与えるだけでなく、血流速度がそれより小さい場合、素直に検出速度ゼロを出力している。このように、有効速度検出範囲において正確な速度検出結果を出力し、かつ、血流速度がそれより小さい場合、素直に、検出速度ゼロを出力するか、あるいは速度検出不能であることを表示することは、本実施形態の速度計測装置の特徴であるということができる。この特徴を用いれば、着目した速度計測装置について、装置内部の構成を仔細に検査するに及ぶことなく、本実施形態の実施の有無を確認することができる。
図３１から図３３には、さらに、クラッタの最終到達速度がナイキスト限界速度の２０％である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。ここまでクラッタ速度が大きくなると、従来方式による速度検出は、完全に破綻している。また、多項式フィット・フィルタ方式は、有効速度検出範囲における誤差が大きくなると同時に、血流速度がそれより小さい場合の誤差も著しく大きくなってしまっている。これに対し、本実施形態の方式は、有効速度検出範囲がやや狭くなるものの、血流速度がナイキスト限界速度の６０％以上の範囲で、ほぼ正確な速度検出結果を与えるだけでなく、血流速度がそれより未満の場合の動作が、やはり素直である。
図３４から図３７には、血流速度自身がナイキスト限界速度の８％に対応するランダムな位相誤差をもつ場合について、それぞれ、式(５)を用いる本実施形態の方式、式(４)を用いる方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。クラッタの最終到達速度については、ナイキスト限界速度の１０％とした。血流速度自身が位相誤差をもつために、負側の低中速領域において、本実施形態の方式も破綻しかかっているが、式(４)を用いる方式では、速度の符号が誤って検出されているのに対し、式(５)を用いる方式では、速度符号については、誤ることなく検出されている。これは、周波数サイドローブレベルを抑圧するよう構築された式(５)を用いる方式の効果である。多項式フィット・フィルタ方式も、同じ負側の低中速領域において、速度検出が破綻しており、その他の速度域においても、かろうじて速度符号が検出できたにとどまっている。また、クラッタ速度が大きいため、従来方式による速度検出は完全に破綻している。
以上述べたように本実施形態によれば、血流を、臓器の運動と峻別して描画することができる。この特徴をさらに生かして、図２３の超音波診断装置では、位相回転検出器１３により検出されたクラッタ信号すなわち臓器エコーの位相回転速度の信号をスキャンコンバータ１９に入力し、臓器の運動速度あるいはその空間微分の分布を示す像を血流像と重畳あるいは並べて表示可能としている。この構成の有用性を、肝腫瘍の場合を例として述べる。肝腫瘍の辺縁部は新生血管が発達し、血流動態が周辺の正常肝とは異なる。また、周辺の正常肝とは固さが異なるため組織の運動も周辺とは異なる。したがって、血流像に加えて臓器の運動速度の場所による変化を表示することは、肝腫瘍の診断にきわめて有用な画像を提供することになる。
また、生体内部の血流など運動する反射物の速度を、そのエコー信号に比べ桁違いに振幅が大きいクラッタ・エコー信号の影響を除いて、正確に検出することができる。具体的には、臓器の超音波探触子への方向の運動速度が、１ｍｍ／ｓ変化するなか、超音波探触子への方向の速度成分 ３ｍｍ／ｓ以上の血流をリアルタイムに描画できる。このように、本実施形態によれば、的確な医療診断を可能とする血流検出・描画機能つき超音波診断装置を提供することが可能となる。すなわち、本実施形態を実施した装置の医用診断上の有用性はきわめて大きく、従って、医用診断を支える工業における本実施形態の意義も、また、大きい。さらに、本実施形態の方法は、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置など、他のパルス・ドップラ装置の運動反射体検出能力をも飛躍的に向上させることができ、本実施形態の工業的ならびに社会的意義は、その意味でもきわめて大である。
（変形例）
本実施形態は前記した実施形態に限定されるものではなく、例えば以下のような種々の変形が可能である。
（１）前記実施形態は、複素展開係数のうち最大の絶対値を持つ係数の次数ｍを求め、次数ｍを境界とした展開係数間の自乗和の比を用いてドップラ速度を算出したが、複素展開係数間の自乗和の比を用いることもできる。
（２）前記実施形態は、超音波を用いて生体内部の血流などの速度分布や速度の空間分布を表示する医療診断用の血流計や血流描画装置に用いたが、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置に用いることができる。
本実施形態は前記した実施形態に限定されるものではなく、例えば以下のような種々の変形が可能である。
（１）前記実施形態は、複素展開係数のうち最大の絶対値を持つ係数の次数ｍを求め、次数ｍを境界とした展開係数間の自乗和の比を用いてドップラ速度を算出したが、複素展開係数間の自乗和の比を用いることもできる。
（２）前記実施形態は、超音波を用いて生体内部の血流などの速度分布や速度の空間分布を表示する医療診断用の血流計や血流描画装置に用いたが、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置に用いることができる。
本発明は、パルス状の反射波を用いた速度計測方法及びこれを用いた速度計測装置に関する。
パルス・ドップラ速度計では、複数のパルス波を運動反射物に送信し、図１のような受信された複数のエコー信号が解析される。すなわち、各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べた時系列信号が解析され、運動反射物に関する速度などの情報が得られる。最も一般的な信号処理方法としては、これらの受信信号を直交検波し、複素時系列信号として解析する。このようなパルス・ドップラ速度計は、生体内部の血流などを検出・描画する超音波診断装置、雨雲などを検出・描画する気象レーダー、飛行物体を検出する航空レーダーなどとして広く用いられている。
もし、反射物が１つしか存在しなければ、前記の時系列信号の位相回転速度 ΔΦ／Δtすなわち符号つき角周波数から、反射物が送受信器に近づく運動速度ｖを、次式のように容易に求めることができる。
ｖ＝λΔΦ／Δt／２ （１）
ここで、λは波長、Δtはパルス送信の時間間隔、ΔΦは位相回転角である。運動速度ｖの符号は、反射物が近づく場合には、運動速度ｖの符号が正となり、遠ざかる場合には負となる。
もし、反射物が１つしか存在しなければ、前記の時系列信号の位相回転速度 ΔΦ／Δtすなわち符号つき角周波数から、反射物が送受信器に近づく運動速度ｖを、次式のように容易に求めることができる。
ｖ＝λΔΦ／Δt／２ （１）
ここで、λは波長、Δtはパルス送信の時間間隔、ΔΦは位相回転角である。運動速度ｖの符号は、反射物が近づく場合には、運動速度ｖの符号が正となり、遠ざかる場合には負となる。
ところで、前記のような実際の受信エコー中では、静止反射物によるエコー信号、いわゆるクラッタ信号の強度が、着目する運動反射物によるエコー信号の強度よりも、一般に、数桁以上大きい。また、静止反射物によるエコー信号は、実際には、時間軸上で完全に静止しているわけではなく、途中の媒質のゆらぎや、医療用超音波診断装置の場合には、血流以外の静止臓器の遅い動きのために、ドリフトしている。このため、静止反射物によるエコー信号の周波数成分は、位相回転速度ゼロの直流成分だけでなく、位相回転速度がゼロでない低周波成分を有している。
このため、前記時系列信号に単純に (１)式の処理を適用しただけでは、運動反射物の速度を検出することはできない。そこで、実際のパルス・ドップラ速度計は、静止反射物によるエコー信号を抑圧し、運動反射物によるエコー信号を相対的に強調するＭＴＩ(Moving Target Indicator) 処理とよばれる信号処理を行った後に、速度の検出あるいは分析を行うように構成されている。
ＭＴＩ処理として最もよく知られているのは、時間領域において畳み込み積和で表現される通常の低域遮断フィルタを用いる処理である。この処理には次の二つの欠点がある。
１） Ｎ１個の時系列データ点を入力とする低域遮断フィルタを用いる場合、後段の速度検出・分析処理部に入力されるデータ点が (Ｎ１−１) 個目減りする。
２） 急峻な遮断特性のフィルタを得にくい。
ＭＴＩ処理として最もよく知られているのは、時間領域において畳み込み積和で表現される通常の低域遮断フィルタを用いる処理である。この処理には次の二つの欠点がある。
１） Ｎ１個の時系列データ点を入力とする低域遮断フィルタを用いる場合、後段の速度検出・分析処理部に入力されるデータ点が (Ｎ１−１) 個目減りする。
２） 急峻な遮断特性のフィルタを得にくい。
パルス・ドップラ速度計の場合は、運動反射物によるエコー信号を保ちながら、静止反射物のドリフト成分を除去する急峻な遮断特性が必要であり、特に、２）が問題となる。また、Ｎ回の送受信により得られるＮ点の時系列信号から速度検出・分析処理をするとき、データ点が目減りする１）の問題は、実際に速度検出・分析演算に用い得るデータ点が(Ｎ−Ｎ１＋１) 点に減ってしまうことを意味する。これは、医療診断用超音波血流描画装置のように実時間性が重要な応用では、望ましいことではない。
この問題を解決する処理としては、非特許文献１に記載されているようにPolynomial Regression Filter （多項式フィット・フィルタ）が提案されている。これは、時系列信号に、０次式（定数）、１次式、２次式、・・・・、Ｍ次式を順次最小２乗フィットし、フィットした成分を引き去ることにより、元の時系列信号のもつドリフト成分を除去する処理である。この処理は、Ｎ点の入力時系列信号にＮ×Ｎの正方行列を掛けることで表現される。このため、出力信号として同じＮ点の時系列信号が得られ、
１） データ点の目減りがない。
また、その低域遮断特性は、フィットする最大次数Ｍの大きさによるが、
２） 同等の遮断周波数をもつ前記の低域遮断フィルタと比較すると、はるかに急峻な遮断特性をもつ。
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Ｖol. 42., 927-937頁 (1995年)
１） データ点の目減りがない。
また、その低域遮断特性は、フィットする最大次数Ｍの大きさによるが、
２） 同等の遮断周波数をもつ前記の低域遮断フィルタと比較すると、はるかに急峻な遮断特性をもつ。
IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Ｖol. 42., 927-937頁 (1995年)
ところが、この多項式フィット・フィルタによる処理の後段において、ｖ＝λΔΦ／Δｔ／２式を用いて運動物体の速度の検出・分析を行うと、遮断周波数の近くで速度演算誤差が、低域遮断フィルタと比較しても大きくなってしまうという問題点がある。このため、多項式フィット・フィルタが有する急峻な遮断特性を十分に使うことができない。
そこで、本発明は、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減することができる速度計測方法及びこれを用いた速度計測装置を提供することを課題とする。
そこで、本発明は、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減することができる速度計測方法及びこれを用いた速度計測装置を提供することを課題とする。
本発明は、速度計測対象物に対してＮ個のパルス波を所定間隔で送信し、受信エコー信号を送信時刻の順番に並べたＮ個の時系列信号を０次から（Ｎ−１）次までの離散的ルジャンドル関数を基底として展開し、この展開係数を用いて速度を算出する。特に、１次以上（Ｎ−１）次以下の（２ｎ−１）次の展開係数と（２ｎ＋１）次の展開係数との線形結合に虚数単位を乗じて、２ｎ次の展開係数と線形結合した２ｎ次の複素展開係数と、（２ｎ＋１）次の展開係数に虚数単位を乗じた後、２ｎ次の展開係数と（２ｎ＋２）次の展開係数と線形結合した（２ｎ＋１）次の複素展開係数とを算出し、各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数ｍを決定し、次数ｍに対応した前記各展開係数間あるいは各複素展開係数間の自乗和の比から運動反射物に関する符号付き速度信号を算出する。
パルス波が照射された運動反射物からは位相回転速度（角周波数）ΔΦ／Δtを有したエコー信号が得られ、ドリフト振動を伴った静止反射物からは低周波成分を有したクラッタ信号が得られる。このクラッタ信号の強度は、エコー信号の強度よりも数桁大きい。エコー信号を次数が大きいほど周波数原点付近で急速にゼロに収束する特性を有する離散的ルジャンドル関数を基底として展開し、各展開係数の横軸を位相回転速度（角周波数）に対応させることによって、運動反射物のエコー信号と静止対象物の低周波成分とが区別される。つまり、静止反射物と運動反射物との反射信号を区別することができる。
本発明によれば、静止反射物と運動反射物との反射信号を区別し、測定誤差を低減できる。
（フィルタの原理）
各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べたＮ点の離散的時系列信号から速度分析を行う場合を考える。このような、パルス・ドップラ速度分析用として典型的な時系列複素信号は、受信されたエコー信号に９０°位相の異なる２つの搬送周波数信号を乗じて得られる１組の信号を、虚数単位を係数として線型結合して複素信号とし、この信号の、各々のパルス送信時刻を基準とする時相の等しい部分を、送信時刻の順番にならべることにより得られる（図１参照）。
このような離散的時系列信号のもつ周波数あるいは位相の分析における常套手段は、離散的フーリエ展開である。具体的には、ｎ＝ １，２，・・・， Ｎ とならんだ時系列信号を、周波数を表す指数ｋ＝−Ｎ／２，・・，０，・・，Ｎ／２の複素正弦波系
Cs(n，k) = cos [ kπ(２ｎ−Ｎ−１)／Ｎ ] + j sin [ kπ(２ｎ−Ｎ−１)/Ｎ ] （２）
で展開する。j は虚数単位である。
負のｋ は反射物が遠ざかる向きの速度に対応し、正のｋ は近づく向きの速度に対応する。ｋ＝±Ｎ／２ はナイキスト限界であり、正負が区別されないので、(２)式により表される (Ｎ＋１)個の複素正弦波関数のうち独立なものはＮ個となり、これらは、直交関数系を形成している。
各々のパルス送信時刻からの経過時間の等しい着目時刻の信号を、送信時刻の順番に並べたＮ点の離散的時系列信号から速度分析を行う場合を考える。このような、パルス・ドップラ速度分析用として典型的な時系列複素信号は、受信されたエコー信号に９０°位相の異なる２つの搬送周波数信号を乗じて得られる１組の信号を、虚数単位を係数として線型結合して複素信号とし、この信号の、各々のパルス送信時刻を基準とする時相の等しい部分を、送信時刻の順番にならべることにより得られる（図１参照）。
このような離散的時系列信号のもつ周波数あるいは位相の分析における常套手段は、離散的フーリエ展開である。具体的には、ｎ＝ １，２，・・・， Ｎ とならんだ時系列信号を、周波数を表す指数ｋ＝−Ｎ／２，・・，０，・・，Ｎ／２の複素正弦波系
Cs(n，k) = cos [ kπ(２ｎ−Ｎ−１)／Ｎ ] + j sin [ kπ(２ｎ−Ｎ−１)/Ｎ ] （２）
で展開する。j は虚数単位である。
負のｋ は反射物が遠ざかる向きの速度に対応し、正のｋ は近づく向きの速度に対応する。ｋ＝±Ｎ／２ はナイキスト限界であり、正負が区別されないので、(２)式により表される (Ｎ＋１)個の複素正弦波関数のうち独立なものはＮ個となり、これらは、直交関数系を形成している。
反射物が１つしか存在せず、それがナイキスト限界範囲内の一定速度で運動している場合を想定し、運動速度の関数として、前記展開により得られる離散的フーリエ展開係数を計算し、Ｎ＝８の場合について、その絶対値を図２にプロットした。横軸は、運動速度に対応する位相回転速度を、パルス送信の時間間隔の逆数 Pulse Repetition Frequency (PRF) の２π倍を単位として表示した。図２では展開係数を表す８個の関数が重なって見づらいので、そのうち４個ずつを見分けやすい線により図３及び図４に示した。
送信時間間隔のＮ倍の時間で、位相が整数回回転するような運動速度では、ゼロでない展開係数はただ１つであり、他の展開係数はすべてゼロとなり、きれいに速度分析されている。ところが、それ以外の位相回転が整数回とならない一般の運動速度では、それぞれの展開係数が、横軸上ピークからはるかに離れた運動速度においても、ピークの１／１０程度の大きさの絶対値をもっている。これは、速度分析計として−２０ｄＢ程度の遠隔的クロストークをもつことを意味する。検出・分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号がドリフトするとき、これはクラッタ信号と運動反射物によるエコー信号とを区別できないという重大な問題を生じさせる。
この問題は、Ｎ個の時系列信号に、ハニング関数のような緩やかに立ち上がって、下がる窓関数による重みをかけることにより、かなり抑圧できるが、本質的には解決しない。また、この方法は、時系列信号の数が大幅に目減りしたのと同様の、望ましくない効果を生じ、特に、時系列信号の数が元々少ない場合には適さない方法である。
具体例として、ナイキスト限界の２/３の一定速度で運動する検出・分析対象である運動反射物があり、その１０００倍のエコー信号振幅をもつクラッタ信号が、ナイキスト限界速度の１／１００でドリフトしている場合を考える。図５には、このときの前記速度分析計の出力である速度スペクトルを示す。検出・分析対象である運動反射物のスペクトルは、速度ゼロ付近にピークをもつクラッタ信号からのクロストーク成分に覆いかくされてしまい、この出力スペクトルの中では検出不可能である。このことが、まさに、クラッタ信号を遮断するＭＴＩフィルタが必要な最も本質的な理由である。
次に、ＭＴＩ処理として最もよく知られ、時間領域において畳み込み積和で表現される通常の低域遮断フィルタについて、その通過振幅特性の典型例を、図６中に実線で示す。このようなフィルタにｎ１ = １,・・・, （Ｎ＋Ｍ−１） とならんだ時系列信号Ｓ１(ｎ１) を、入力すると、出力としてｎ０＝１,・・・, Ｎ とならんだ時系列信号Ｓ０(ｎ０)が得られる。
Ｓ０(ｎ０)＝ΣＦ（ｍ）Ｓ１(ｎ０＋ｍ−１) (３)
ここで、Σは ｍ＝１，・・・，Ｍの和を表す。図６では、低域遮断フィルタとして最も簡単なＭ＝３，Ｆ(１)＝−１， Ｆ（２）＝２, Ｆ（３）＝−１の場合について実線で示した。図５と同様の入力信号を、この低域遮断フィルタを通した後、図５と同様に速度分析すると、図７中の実線のようになる。低域遮断フィルタの働きによりクラッタ信号の振幅が１／２０００倍程度に抑圧された結果、ナイキスト限界の０．３倍から０．４倍の速度にピークをもつ運動反射物のスペクトルが見えるようになっている。
Ｓ０(ｎ０)＝ΣＦ（ｍ）Ｓ１(ｎ０＋ｍ−１) (３)
ここで、Σは ｍ＝１，・・・，Ｍの和を表す。図６では、低域遮断フィルタとして最も簡単なＭ＝３，Ｆ(１)＝−１， Ｆ（２）＝２, Ｆ（３）＝−１の場合について実線で示した。図５と同様の入力信号を、この低域遮断フィルタを通した後、図５と同様に速度分析すると、図７中の実線のようになる。低域遮断フィルタの働きによりクラッタ信号の振幅が１／２０００倍程度に抑圧された結果、ナイキスト限界の０．３倍から０．４倍の速度にピークをもつ運動反射物のスペクトルが見えるようになっている。
同じ目的のＭＴＩフィルタとして、前記の多項式フィット・フィルタを設計したときの通過振幅特性の典型例を、図６中に破線で示す。この例では、０次から３次式までをフィットして引き去った。前記の畳み込み型フィルタよりも急峻な低域遮断特性が見てとれる。また、図５と同様の入力信号を、この多項式フィット・フィルタを通した後、図５と同様に速度分析した結果を、図７中に破線で示した。低域遮断特性の急峻さのために、クラッタ信号が完全に抑圧されている。多項式フィット・フィルタ特有の急峻な低域遮断特性は以下のように理解することができる。
連続座標系における狭義ルジャンドル多項式については、岩波数学公式III特殊関数８２−８５頁 に記述されており、例えば、８次まで記すと
Ｐ0（ｘ）＝１
Ｐ1（ｘ）＝ｘ＝cosθ
Ｐ2（ｘ）＝（１／２）（３ｘ2−１）＝（１／４）（３cos２θ＋１）
Ｐ3（ｘ）＝（１／２）（５ｘ3−３ｘ）＝（１／８）（５cos３θ＋３cosθ）
Ｐ4（ｘ）＝（１／８）（３５ｘ4−３０ｘ2＋３）＝（１／６４）（３５cos４θ＋２０cos２θ＋９）
Ｐ5（ｘ）＝（１／８）（６３ｘ5−７０ｘ3＋１５ｘ）＝（１／１２８）（６３cos５θ＋３５cos３θ＋３０cosθ）
Ｐ6（ｘ）＝（１／１６）（２３１ｘ6−３１５ｘ4＋１０５ｘ2−５）＝（１／５１２）（２３１cos６θ＋１２６cos４θ＋１０５cos２θ＋５０）
Ｐ7（ｘ）＝（１／１６）（４２９ｘ7−６３９ｘ5＋３１５ｘ3−３５ｘ）＝（１／２10）（４２９cos７θ＋２３１cos５θ＋１８９cos３θ＋１７５cosθ）
Ｐ8（ｘ）＝（１／１２８）（６４３５ｘ8−１２０１２ｘ6＋６９３０ｘ4−１２６０ｘ2＋３５）＝（１／２14）（６４３５cos８θ＋３４３２cos６θ＋２７７２cos４θ＋２５２０cos２θ＋１２２５）
である。離散的座標系についても使えるように、一般化すると、偶数（２ｎ）次ルジャンドル多項式は、２ｎより低い次数のルジャンドル多項式すべてと直交する最高次数２ｎの偶関数と定義することができ、奇数（２ｎ＋１）次ルジャンドル多項式は、(２ｎ＋１)より低い次数のルジャンドル多項式すべてと直交する最高次数(２ｎ＋１)の奇関数と定義することができる。
Ｐ0（ｘ）＝１
Ｐ1（ｘ）＝ｘ＝cosθ
Ｐ2（ｘ）＝（１／２）（３ｘ2−１）＝（１／４）（３cos２θ＋１）
Ｐ3（ｘ）＝（１／２）（５ｘ3−３ｘ）＝（１／８）（５cos３θ＋３cosθ）
Ｐ4（ｘ）＝（１／８）（３５ｘ4−３０ｘ2＋３）＝（１／６４）（３５cos４θ＋２０cos２θ＋９）
Ｐ5（ｘ）＝（１／８）（６３ｘ5−７０ｘ3＋１５ｘ）＝（１／１２８）（６３cos５θ＋３５cos３θ＋３０cosθ）
Ｐ6（ｘ）＝（１／１６）（２３１ｘ6−３１５ｘ4＋１０５ｘ2−５）＝（１／５１２）（２３１cos６θ＋１２６cos４θ＋１０５cos２θ＋５０）
Ｐ7（ｘ）＝（１／１６）（４２９ｘ7−６３９ｘ5＋３１５ｘ3−３５ｘ）＝（１／２10）（４２９cos７θ＋２３１cos５θ＋１８９cos３θ＋１７５cosθ）
Ｐ8（ｘ）＝（１／１２８）（６４３５ｘ8−１２０１２ｘ6＋６９３０ｘ4−１２６０ｘ2＋３５）＝（１／２14）（６４３５cos８θ＋３４３２cos６θ＋２７７２cos４θ＋２５２０cos２θ＋１２２５）
である。離散的座標系についても使えるように、一般化すると、偶数（２ｎ）次ルジャンドル多項式は、２ｎより低い次数のルジャンドル多項式すべてと直交する最高次数２ｎの偶関数と定義することができ、奇数（２ｎ＋１）次ルジャンドル多項式は、(２ｎ＋１)より低い次数のルジャンドル多項式すべてと直交する最高次数(２ｎ＋１)の奇関数と定義することができる。
例として、一連の８点からなる時系列信号について、０次から７次までの離散的ルジャンドル関数を計算し、偶数次を図８に、奇数次を図９に示した。０次・１次を実線、２次・３次を破線、４次・５次を１点鎖線、６次・７次を点線にて示した。振幅は、ＲＭＳ（２乗平均の平方根）値が１となるよう規格化した。符号は、最低次項の符号が正となるよう、すなわち、偶数次多項式については定数項、奇数次多項式については１次項の符号が正となるよう定めた。図１０には、その周波数スペクトルをプロットした。周波数成分の振幅の絶対値を、個々のルジャンドル展開係数の最大値で規格化して示した。図１０では８個の関数が重なってしまって見づらいので、そのうち４個ずつを見分けやすい線により図１１及び図１２に示した。Ａ（０），Ａ（７）を実線、Ａ（１），Ａ（６）を破線、Ａ（２），Ａ（５）を1点鎖線、Ａ（３），Ａ（４）を点線にて示した。
図１０から図１２を、図２から図４と同様、入力複素正弦波に対する周波数応答として眺めてみる。ピークを主応答、それ以外を副応答とみる。ルジャンドル関数は、図２から図４のフーリエ展開系に比べ、主応答より高周波側の副応答は大きいが、低周波側には副応答をもたず、次数が大きいほど周波数原点付近で急速にゼロに収束する。ここでは解析的証明は省略するが、周波数の関数としてのｍ次ルジャンドル関数は、周波数原点において０次から(ｍ−１)次導関数まですべてが０となる。
前記の多項式フィット・フィルタは、入力信号を、この離散的ルジャンドル関数を基底として展開し、０からＭ次までを除き、(Ｍ＋１)次以上の成分の和を出力信号とするフィルタであるということができる。したがって、出力信号は、周波数原点付近で周波数の(Ｍ＋１)乗に比例する挙動でゼロに収束していく。図６及び図７中に破線で示された例のように、多項式フィット・フィルタの低域遮断特性が優れているのは、このようなルジャンドル関数の本質的な特性に由来するものである。このように考えると、クラッタという周波数原点付近に桁違いに大きな振幅をもつ妨害成分の影響を抑えて、着目する運動物体を検出したり、その速度を分析したりしなければならないドップラ速度検出には、多項式フィット・フィルタやルジャンドル関数展開を利用する方式は、本質的に適しているといえる。
フーリエ展開系では、複素正弦波を基底とすることにより、図２から図５や図７に示したように符号つきの速度分析が可能であることは広く知られており、そのような分析方法が広く用いられている。ところが、ルジャンドル展開系は、前記のように、ドップラ速度検出に極めて適した性質をもつものの、このまま用いたのでは、図１０から図２０より見てとれるように、絶対値が等しいが符号の異なる速度に対して等しい応答を示すので、符号つきの速度分析は不可能である。
そこで、次数が１つだけ異なる偶数次と奇数次のルジャンドル関数を虚数単位 j を係数として線型結合して、複素ルジャンドル関数を形成し、これにより符号つきの速度分析を可能とする。これは、偶数次のルジャンドル関数が余弦波に、奇数次のルジャンドル関数が正弦波に似た変化を示すことに着目し、余弦関数と正弦関数を虚数単位 j を係数として線型結合することにより、位相角をｊ倍した指数関数、すなわち複素正弦関数が得られ、この位相角の増減から符号つき速度が得られるのと同様に、複素ルジャンドル関数から符号つき速度を得ようとするものである。
具体的には、一連のＮ点からなる時系列信号を、まず、０次から (Ｎ−１)次の離散的ルジャンドル関数を基底として展開して、Ａ（０）からＡ（Ｎ−１）までの展開係数を得て、これをもとに、一連の複素ルジャンドル係数
Ｃ（±(２ｎ−３／２)）＝Ａ（２ｎ−２）±ｊＡ（２ｎ−１） (１≦２ｎ−１≦Ｎ-1)
Ｃ（±(２ｎ−１／２)）＝Ａ（２ｎ）±j Ａ（２ｎ−１） (２≦２ｎ≦Ｎ−1;ｎ:自然数, 複号同順) (４)
を算出し、この係数の相対的な大きさにより符号つき速度分析を行う。
具体的には、一連のＮ点からなる時系列信号を、まず、０次から (Ｎ−１)次の離散的ルジャンドル関数を基底として展開して、Ａ（０）からＡ（Ｎ−１）までの展開係数を得て、これをもとに、一連の複素ルジャンドル係数
Ｃ（±(２ｎ−３／２)）＝Ａ（２ｎ−２）±ｊＡ（２ｎ−１） (１≦２ｎ−１≦Ｎ-1)
Ｃ（±(２ｎ−１／２)）＝Ａ（２ｎ）±j Ａ（２ｎ−１） (２≦２ｎ≦Ｎ−1;ｎ:自然数, 複号同順) (４)
を算出し、この係数の相対的な大きさにより符号つき速度分析を行う。
図１３から図１５には、図１０から図１２と同様Ｎ＝８の場合について、複素正弦波入力に対する出力としての複素ルジャンドル係数の絶対値を、符号つき速度に対応する位相回転速度の関数として示した。複素ルジャンドル係数の絶対値は、それぞれの最大値で規格化してプロットした。図１３中の実線はＣ（６．５）、破線はＣ（−６．５）、図１４中の実線はＣ（５．５）、破線はＣ（−５．５）、図１５中の実線はＣ（４．５）、破線はＣ（−４．５）である。これらの図から、複素ルジャンドル係数が、速度に対し正負を区別して応答していることがわかる。なお、このような速度に対する応答をもつ複素ルジャンドル系がＮ≦３５の範囲で構築できることが確認されている。
図１３から図１５に示した複素ルジャンドル係数の応答は、ルジャンドル展開係数の特徴である位相回転速度原点すなわちドップラ周波数原点のまわりの極めて優れたカットオフ特性を示している反面、図２から図４に示したフーリエ展開係数の応答と比較して、周波数サイドローブレベルが高めである。すなわち、周波数サイドローブレベルは、フーリエ展開係数の図２から図４の例では、０．２程度であるのに対し、複素ルジャンドル係数の図１４の例では、０．６近くに達している。周波数サイドローブレベルが高いことは、Ｓ/Ｎ比すなわちノイズレベルに対する信号レベルが低い場合に、速度検出誤差を生む可能性が高いことを意味するので、周波数サイドローブレベルは抑圧することが望ましい。
そこで、Ａ（２ｎ）をもとにＣ（±２ｎ）を生成する際、Ａ（２ｎ）と同じ周波数において最大値をとるようＡ（２ｎ−１）とＡ（２ｎ＋１）とを線型結合して用い、また、Ａ（２ｎ−１）をもとにＣ（±(２ｎ−１)）を生成する際、Ａ（２ｎ−１）と同じ周波数において最大値をとるようＡ（２ｎ−２）とＡ（２ｎ）とを線型結合して用いる。
すなわち、
Ｃ（±(２ｎ−１)）＝α(２ｎ−１)×Ａ（２ｎ−２) ＋β(２ｎ−１)×Ａ（２ｎ）±j Ａ（２ｎ−１） (１≦２ｎ−１≦Ｎ−１)
Ｃ（±２ｎ）＝Ａ（２ｎ）±j [α(２ｎ)×Ａ（２ｎ−１）+ β(２ｎ)×Ａ（２ｎ＋１） ] (２≦２ｎ≦Ｎ−１；ｎ：自然数，複号同順；α及びβ：正の実数 ) (５)
となる。
すなわち、
Ｃ（±(２ｎ−１)）＝α(２ｎ−１)×Ａ（２ｎ−２) ＋β(２ｎ−１)×Ａ（２ｎ）±j Ａ（２ｎ−１） (１≦２ｎ−１≦Ｎ−１)
Ｃ（±２ｎ）＝Ａ（２ｎ）±j [α(２ｎ)×Ａ（２ｎ−１）+ β(２ｎ)×Ａ（２ｎ＋１） ] (２≦２ｎ≦Ｎ−１；ｎ：自然数，複号同順；α及びβ：正の実数 ) (５)
となる。
図１６から図１７には、図１３から図１５と同様Ｎ＝８の場合について、ルジャンドル係数Ａ（４）、Ａ（５）、Ａ（６）から複素ルジャンドル係数Ｃ（±５）を生成する例を示した。図１６中の実線はＡ（５）の絶対値のドップラ周波数応答である。一方、点線は、そのＡ（５）と同じ周波数において最大値をとるようにＡ（４）とＡ（６）とを線型結合して生成した係数の絶対値のドップラ周波数応答である。図１７は、そのＡ（４）とＡ（６）との線型結合を実部、Ａ（５）を虚部として得られる複素ルジャンドル係数Ｃ（±５）の絶対値のドップラ周波数応答である。図中実線はＣ（５）であり、破線はＣ（−５）である。周波数サイドローブレベルは、０．４程度に抑えられている。これは、図１４に示されたＣ（±５．５）や図１５に示されたＣ（±４．５）の周波数サイドローブレベル０．６程度に比べ、大きく改善されている。同様に生成した複素ルジャンドル係数Ｃ（±６）、Ｃ（±４）、Ｃ（±３）の絶対値のドップラ周波数応答を図１８から図２０に示した。これらは図１７と同様に、実線がＣ（６）、Ｃ（４）、Ｃ（３）であり、破線がＣ（−６）、Ｃ（−４）、Ｃ（−３）である。周波数サイドローブレベルは、やはり０．４程度に抑えられている。すなわち、Ａ（２ｎ−１）と同じ角周波数において最大値をとるようＡ（２ｎ−２）とＡ（２ｎ）とを線型結合することによって、サイドローブレベルが０．６から０．４程度に抑えられる。
図２１には、ｍを自然数とするとき、(ｍ−１)次から(Ｎ−１)次までのルジャンドル係数の自乗和とｍ次から(Ｎ−１)次までのルジャンドル係数の自乗和との比
Ｒ（ｍ）＝[Ａ（ｍ）2＋・・・＋Ａ（Ｎ−１)2] ／ [Ａ（ｍ−１）2＋Ａ（ｍ）2＋・・・＋Ａ（Ｎ−１)2] (１≦ｍ≦Ｎ−１) (６)
を、前記の例と同様Ｎ＝８の場合について示した。これらの曲線は、ドップラ周波数原点から始まり、原点から離れるにともない単調増加して、一旦最大値をとり、その後リップルする。その単調増加部分を用いることにより、図２２に示すような、速度校正曲線を得ることができる。絶対値が最大の複素ルジャンドル係数の次数が求まると、その次数に対応する速度校正曲線を用いて、ルジャンドル係数の自乗和の比からドップラ速度を算出すればよい。
Ｒ（ｍ）＝[Ａ（ｍ）2＋・・・＋Ａ（Ｎ−１)2] ／ [Ａ（ｍ−１）2＋Ａ（ｍ）2＋・・・＋Ａ（Ｎ−１)2] (１≦ｍ≦Ｎ−１) (６)
を、前記の例と同様Ｎ＝８の場合について示した。これらの曲線は、ドップラ周波数原点から始まり、原点から離れるにともない単調増加して、一旦最大値をとり、その後リップルする。その単調増加部分を用いることにより、図２２に示すような、速度校正曲線を得ることができる。絶対値が最大の複素ルジャンドル係数の次数が求まると、その次数に対応する速度校正曲線を用いて、ルジャンドル係数の自乗和の比からドップラ速度を算出すればよい。
検出・分析対象である運動反射物エコーに比べて桁違いに大きな振幅をもつクラッタ信号が存在し、それがドリフトするときには、そのドリフトの程度が大きくなるに応じて、Ｃ（±１）、Ｃ（±２）、……と、次数の小さな係数を順次除いた複素ルジャンドル係数系を用いて、速度検出・分析を行えばよい。このようにルジャンドル係数の遮断次数を制御することにより、クラッタ信号の影響を抑圧しながら、運動反射物の速度検出・分析を行うことができる。
（第１実施形態）
以下、本発明の実施形態について説明する。
図２３は、本発明の一実施形態である速度計測装置を使用した超音波診断装置のブロック図であり、血流描画機能を備えている。超音波探触子１を構成する各素子は、切り替えスイッチ群２を介して、送波ビームフォーマ３と受波ビームフォーマ１０に接続されている。送波ビームフォーマ３では、送受信シークエンス制御部６による制御に従って、送信波形メモリ５から送信波形選択部４により選択されて読み出された波形を用い、各素子を通じて送信されたときに指向性を持つ超音波パルスとなるような信号が生成される。この信号が、超音波探触子１の各素子により超音波パルスに変換されて生体に送信される。生体中で反射あるいは散乱されて超音波探触子１に戻ってきた超音波エコー信号は、各素子に受信されて、電気信号に変換される。受波ビームフォーマ１０では、送受信シークエンス制御部６による制御に従って、指向性を持つ受信感度を生成すべく、各受波信号に遅延時間を与えて互いに加算する。遅延加算により得られた時系列信号は、やはり送受信シークエンス制御部６による制御に従って、受波メモリ選択部１１により選択された受波信号メモリ１２中の１つのバンクへ一旦書き込まれ、ドップラ信号分析すべきＮ個の時系列信号がそろったのちに読み出されて、位相回転検出器１３，１５、ダウン・ミキシング器１４、血流信号検出・分析部１６において速度を検出・分析すべく信号処理される。
以下、本発明の実施形態について説明する。
図２３は、本発明の一実施形態である速度計測装置を使用した超音波診断装置のブロック図であり、血流描画機能を備えている。超音波探触子１を構成する各素子は、切り替えスイッチ群２を介して、送波ビームフォーマ３と受波ビームフォーマ１０に接続されている。送波ビームフォーマ３では、送受信シークエンス制御部６による制御に従って、送信波形メモリ５から送信波形選択部４により選択されて読み出された波形を用い、各素子を通じて送信されたときに指向性を持つ超音波パルスとなるような信号が生成される。この信号が、超音波探触子１の各素子により超音波パルスに変換されて生体に送信される。生体中で反射あるいは散乱されて超音波探触子１に戻ってきた超音波エコー信号は、各素子に受信されて、電気信号に変換される。受波ビームフォーマ１０では、送受信シークエンス制御部６による制御に従って、指向性を持つ受信感度を生成すべく、各受波信号に遅延時間を与えて互いに加算する。遅延加算により得られた時系列信号は、やはり送受信シークエンス制御部６による制御に従って、受波メモリ選択部１１により選択された受波信号メモリ１２中の１つのバンクへ一旦書き込まれ、ドップラ信号分析すべきＮ個の時系列信号がそろったのちに読み出されて、位相回転検出器１３，１５、ダウン・ミキシング器１４、血流信号検出・分析部１６において速度を検出・分析すべく信号処理される。
読み出されたＮ個の時系列受信信号は、各々を得るために行われたパルス送信の時刻を基準とする時相の等しい部分を、送信時刻の順番に並べたのちに処理される。この速度検出アルゴリズムについて、図２４のフローチャート及び図２３を参照して説明する。まず、位相回転検出器１３では、クラッタ信号を含んだままの時系列信号から位相回転を求める。位相回転を検出する最も典型的な方法はｎ＝１,・・・, Ｎとならんだ隣接する時系列複素信号 Ｓ（ｎ），Ｓ（ｎ＋１）をもとに、
Ｐ（ｎ）＝Ｓ（ｎ＋１）Ｓ（ｎ）* ／||Ｓ（ｎ＋１）||／||Ｓ（ｎ）|| (ｎ＝１,・・・,Ｎ−１) (７)
を計算し、Ｐ（ｎ）の平均値Ｐａの位相から平均位相回転速度を求める。
ここで、Ｓ（ｎ）*はＳ（ｎ）の複素共役、||Ｓ（ｎ）|| はＳ（ｎ）の絶対値である。一般に、クラッタ・エコー信号は血流エコー信号より桁違いに信号振幅が大きいので、Ｐａは、概ねクラッタ・エコー信号の位相回転平均値と考えてよい。すなわち、隣接サンプル間位相差平均値よりクラッタ平均速度が計算される（図２４のＳ１）。また、時系列複素信号Ｓ（ｎ）は、受信されたエコー信号に９０°位相の異なる２つの搬送周波数信号を乗じて得られる。
Ｐ（ｎ）＝Ｓ（ｎ＋１）Ｓ（ｎ）* ／||Ｓ（ｎ＋１）||／||Ｓ（ｎ）|| (ｎ＝１,・・・,Ｎ−１) (７)
を計算し、Ｐ（ｎ）の平均値Ｐａの位相から平均位相回転速度を求める。
ここで、Ｓ（ｎ）*はＳ（ｎ）の複素共役、||Ｓ（ｎ）|| はＳ（ｎ）の絶対値である。一般に、クラッタ・エコー信号は血流エコー信号より桁違いに信号振幅が大きいので、Ｐａは、概ねクラッタ・エコー信号の位相回転平均値と考えてよい。すなわち、隣接サンプル間位相差平均値よりクラッタ平均速度が計算される（図２４のＳ１）。また、時系列複素信号Ｓ（ｎ）は、受信されたエコー信号に９０°位相の異なる２つの搬送周波数信号を乗じて得られる。
ダウン・ミキシング器１４では、ここで求めた平均値Ｐａ又はＰａの周辺領域における空間平均値を用い、エコー信号の位相回転平均値がゼロとなるようにミキシング処理が行われる。いいかえれば、ダウン・ミキシングによりクラッタ平均速度による位相回転が打ち消される（図２４のＳ２参照）。すなわち、Ｓ（ｎ）とＰａより
Ｓｄ（ｎ）＝Ｓ（ｎ）Ｐａ* n (ｎ＝１,・・・, Ｎ) (８)
を得る。このダウン・ミキシング処理を行うことにより、クラッタ信号を抑圧する後段の処理をより効果的に行うことができる。
Ｓｄ（ｎ）＝Ｓ（ｎ）Ｐａ* n (ｎ＝１,・・・, Ｎ) (８)
を得る。このダウン・ミキシング処理を行うことにより、クラッタ信号を抑圧する後段の処理をより効果的に行うことができる。
得られた隣接サンプルである時系列複素信号Ｓｄ（ｎ），Ｓｄ（ｎ＋１）をもとに、位相回転検出器１５では、位相回転検出器１３と同様に、
Ｐｄ（ｎ）＝Ｓｄ（ｎ＋１）Ｓｄ（ｎ）* ／ ||Ｓｄ（ｎ＋１）|| ／||Ｓｄ（ｎ）|| (ｎ＝１,・・・, Ｎ−１) (９)
を計算する。この隣接サンプル間位相差最大値より、消し残したクラッタ速度の最大値が評価される（図２４のＳ３）。この位相回転最大値の大きさに応じて、血流信号検出・分析部１６の遮断特性を制御するために、ルジャンドル展開係数の最低次数が決定される（図２４のＳ４）。すなわち、位相回転最大値が大きいときには、小さいときに比べ、ルジャンドル展開係数の遮断次数Ｍがより高く設定され、クラッタ成分が効果的に抑圧される。
Ｐｄ（ｎ）＝Ｓｄ（ｎ＋１）Ｓｄ（ｎ）* ／ ||Ｓｄ（ｎ＋１）|| ／||Ｓｄ（ｎ）|| (ｎ＝１,・・・, Ｎ−１) (９)
を計算する。この隣接サンプル間位相差最大値より、消し残したクラッタ速度の最大値が評価される（図２４のＳ３）。この位相回転最大値の大きさに応じて、血流信号検出・分析部１６の遮断特性を制御するために、ルジャンドル展開係数の最低次数が決定される（図２４のＳ４）。すなわち、位相回転最大値が大きいときには、小さいときに比べ、ルジャンドル展開係数の遮断次数Ｍがより高く設定され、クラッタ成分が効果的に抑圧される。
ここで、本実施形態の特徴部である血流信号検出・分析部１６の動作について、もう少し詳しく述べる。入力された時系列信号 Ｓｄ(１) ,・・・, Ｓｄ（Ｎ） は、まず、図８及び図９のような０から(Ｎ−１)次の離散的ルジャンドル関数に展開され、その展開係数Ａ(０) ,・・・, Ａ(Ｎ−１) が求められる。この計算は、次のような行列演算により容易に行うことができる。ｎ次ルジャンドル関数を行ベクトル Ｌ(ｎ) で表し、これをｎ＝０から(Ｎ−１)までならべてＮ×Ｎ行列ＬＬを生成すると、Ｓｄ(１) ,・・・, Ｓｄ(Ｎ)を要素とする列ベクトルＳｄからＡ(０) ,・・・,Ａ(Ｎ−１)を要素とする列ベクトルＡを得るための行列は、
Ｆ＝(Ｌ・tＬＬ)-1ＬＬ （１０）
と求めることができる。ここで、tＬＬはＬＬの転置行列、ＬＬ-1はＬＬの逆行列を表す。この行列を予め準備しておけば、行列演算
Ａ＝Ｆ・Ｓｄ (１１)
によって迅速に、ルジャンドル展開係数を求めることができる。この展開係数Ａ（０）,・・・, Ａ（Ｎ−１）から、式(５)に従って、Ｃ（±１) ,・・・, Ｃ（±(Ｎ−１))の複素ルジャンドル係数が生成される（図２４のＳ５）。いいかえれば、ｎ次ルジャンドル関数を行ベクトルＬ（ｎ）で表し、これをｎ＝０から（Ｎ−１）まで列方向に並べたＮ×Ｎ行列ＬＬを転置した転置行列 ｔ ＬＬを生成する過程と、転置行列 ｔ ＬＬにＮ×Ｎ行列ＬＬを前から乗算したものの逆行列（ＬＬ・ ｔ ＬＬ)-1を算出する過程と、Ｎ×Ｎ行列ＬＬに逆行列（ＬＬ・ ｔ ＬＬ)-1を前から乗算する過程とによって算出された行列ＦＦを要素が時系列信号である列ベクトルＳｄ（ｎ）に前から乗算して、ルジャンドル展開係数を求めることができる。前記したように別途、位相回転最大値に応じてルジャンドル展開係数の最低次数（遮断次数Ｍ）が決定され（図２４のＳ４）、Ｃ（±１) から Ｃ（±(Ｍ＋１)) までの低次数複素ルジャンドル係数が棄却される。
Ｆ＝(Ｌ・tＬＬ)-1ＬＬ （１０）
と求めることができる。ここで、tＬＬはＬＬの転置行列、ＬＬ-1はＬＬの逆行列を表す。この行列を予め準備しておけば、行列演算
Ａ＝Ｆ・Ｓｄ (１１)
によって迅速に、ルジャンドル展開係数を求めることができる。この展開係数Ａ（０）,・・・, Ａ（Ｎ−１）から、式(５)に従って、Ｃ（±１) ,・・・, Ｃ（±(Ｎ−１))の複素ルジャンドル係数が生成される（図２４のＳ５）。いいかえれば、ｎ次ルジャンドル関数を行ベクトルＬ（ｎ）で表し、これをｎ＝０から（Ｎ−１）まで列方向に並べたＮ×Ｎ行列ＬＬを転置した転置行列 ｔ ＬＬを生成する過程と、転置行列 ｔ ＬＬにＮ×Ｎ行列ＬＬを前から乗算したものの逆行列（ＬＬ・ ｔ ＬＬ)-1を算出する過程と、Ｎ×Ｎ行列ＬＬに逆行列（ＬＬ・ ｔ ＬＬ)-1を前から乗算する過程とによって算出された行列ＦＦを要素が時系列信号である列ベクトルＳｄ（ｎ）に前から乗算して、ルジャンドル展開係数を求めることができる。前記したように別途、位相回転最大値に応じてルジャンドル展開係数の最低次数（遮断次数Ｍ）が決定され（図２４のＳ４）、Ｃ（±１) から Ｃ（±(Ｍ＋１)) までの低次数複素ルジャンドル係数が棄却される。
次に、棄却されないで残された正次数の複素ルジャンドル係数の自乗和と負次数の複素ルジャンドル係数の自乗和とを比較して、速度符号が決定される（図２４のＳ６参照）。具体的に、大きな方の次数の符号を速度の符号とする。さらに、速度符号次数の複素ルジャンドル係数のうち、絶対値最大の係数が検出される（図２４のＳ７参照）。具体的に、正なら正の次数、負なら負の次数の複素ルジャンドル係数と、最大次数(Ｎ−１)のルジャンドル係数のうち、絶対値最大の係数の次数ｍが検出される。
そして、検出したルジャンドル係数に対応する校正曲線を用いて血流速度を算出する（図２４のＳ８参照）。すなわち、次数ｍに対応する図２２に示したような速度校正曲線を用いて、式(６)のように(ｍ−１)次から(Ｎ−１)次までのルジャンドル係数の自乗和とｍ次から(Ｎ−１)次までのルジャンドル係数の自乗和との比からドップラ速度を求める。
得られた血流速度信号は、エコー振幅検出器１７及びエコー振幅圧縮器１８により得られる静止臓器からのエコー信号と共に、スキャンコンバータ１９に入力される。スキャンコンバータでは、入力された複数の信号を適宜重畳して表示器２０にて２次元ないし３次元表示すべく、信号の生成・制御を行う。
次に、本実施形態のドップラ速度計測装置の動作例を従来方式と比較して示す。クラッタ・エコー信号の振幅は、血流エコーの３００倍とし、クラッタ速度については、初期値がゼロで一定の加速度で立ち上がる場合を想定した。時系列信号数Ｎ＝８の場合について、数値計算シミュレーションを行って比較した。従来方式としては、ＭＴＩフィルタとして図６の場合の低域遮断フィルタを用い、その後段において、式(５)を用いて位相回転量平均値Ｐａを得たのと同様の方法で、速度を算出する処理を血流信号検出・分析部１６において行う例を選んだ。さらに、ＭＴＩフィルタを、やはり図５に例を示したような多項式フィット・フィルタに置き換えた場合とも、比較した。この場合、遮断次数Mは、本発明の方法と同様、位相回転検出器１５の出力信号に応じて制御した。
図２５から図２７には、クラッタの最終到達速度がナイキスト限界速度の０．８％である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、及び多項式フィット・フィルタ方式による速度検出結果を示した。横軸に入力としての血流速度、縦軸に出力として検出された速度をとり、有効速度検出範囲を実線、範囲外を点線によりプロットした。また、理想的場合を一点鎖線により示した。このように、クラッタ速度が低い場合には、いずれの方式も、ほぼ正確な速度分析結果を与えているが、本実施形態の方式による誤差が最も小さい。
図２８から図３０には、同様に、クラッタの最終到達速度がナイキスト限界速度の３％である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。このクラッタ速度域では、低域遮断フィルタ方式による速度検出は、ほぼ破綻に至っている。また、多項式フィット・フィルタ方式は、有効速度検出範囲における誤差は大きくないものの、血流速度がそれより小さい場合の誤差が大きく、何らかの対策を加えない限り、使うことができない。これに対し、本実施形態の方式は、有効速度検出範囲において正確な速度検出結果を与えるだけでなく、血流速度がそれより小さい場合、素直に検出速度ゼロを出力している。このように、有効速度検出範囲において正確な速度検出結果を出力し、かつ、血流速度がそれより小さい場合、素直に、検出速度ゼロを出力するか、あるいは速度検出不能であることを表示することは、本実施形態の速度計測装置の特徴であるということができる。この特徴を用いれば、着目した速度計測装置について、装置内部の構成を仔細に検査するに及ぶことなく、本実施形態の実施の有無を確認することができる。
図３１から図３３には、さらに、クラッタの最終到達速度がナイキスト限界速度の２０％である場合について、それぞれ、本実施形態の方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。ここまでクラッタ速度が大きくなると、従来方式による速度検出は、完全に破綻している。また、多項式フィット・フィルタ方式は、有効速度検出範囲における誤差が大きくなると同時に、血流速度がそれより小さい場合の誤差も著しく大きくなってしまっている。これに対し、本実施形態の方式は、有効速度検出範囲がやや狭くなるものの、血流速度がナイキスト限界速度の６０％以上の範囲で、ほぼ正確な速度検出結果を与えるだけでなく、血流速度がそれより未満の場合の動作が、やはり素直である。
図３４から図３７には、血流速度自身がナイキスト限界速度の８％に対応するランダムな位相誤差をもつ場合について、それぞれ、式(５)を用いる本実施形態の方式、式(４)を用いる方式、低域遮断フィルタ方式、多項式フィット・フィルタ方式による速度検出結果を示した。クラッタの最終到達速度については、ナイキスト限界速度の１０％とした。血流速度自身が位相誤差をもつために、負側の低中速領域において、本実施形態の方式も破綻しかかっているが、式(４)を用いる方式では、速度の符号が誤って検出されているのに対し、式(５)を用いる方式では、速度符号については、誤ることなく検出されている。これは、周波数サイドローブレベルを抑圧するよう構築された式(５)を用いる方式の効果である。多項式フィット・フィルタ方式も、同じ負側の低中速領域において、速度検出が破綻しており、その他の速度域においても、かろうじて速度符号が検出できたにとどまっている。また、クラッタ速度が大きいため、従来方式による速度検出は完全に破綻している。
以上述べたように本実施形態によれば、血流を、臓器の運動と峻別して描画することができる。この特徴をさらに生かして、図２３の超音波診断装置では、位相回転検出器１３により検出されたクラッタ信号すなわち臓器エコーの位相回転速度の信号をスキャンコンバータ１９に入力し、臓器の運動速度あるいはその空間微分の分布を示す像を血流像と重畳あるいは並べて表示可能としている。この構成の有用性を、肝腫瘍の場合を例として述べる。肝腫瘍の辺縁部は新生血管が発達し、血流動態が周辺の正常肝とは異なる。また、周辺の正常肝とは固さが異なるため組織の運動も周辺とは異なる。したがって、血流像に加えて臓器の運動速度の場所による変化を表示することは、肝腫瘍の診断にきわめて有用な画像を提供することになる。
また、生体内部の血流など運動する反射物の速度を、そのエコー信号に比べ桁違いに振幅が大きいクラッタ・エコー信号の影響を除いて、正確に検出することができる。具体的には、臓器の超音波探触子への方向の運動速度が、１ｍｍ／ｓ変化するなか、超音波探触子への方向の速度成分 ３ｍｍ／ｓ以上の血流をリアルタイムに描画できる。このように、本実施形態によれば、的確な医療診断を可能とする血流検出・描画機能つき超音波診断装置を提供することが可能となる。すなわち、本実施形態を実施した装置の医用診断上の有用性はきわめて大きく、従って、医用診断を支える工業における本実施形態の意義も、また、大きい。さらに、本実施形態の方法は、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置など、他のパルス・ドップラ装置の運動反射体検出能力をも飛躍的に向上させることができ、本実施形態の工業的ならびに社会的意義は、その意味でもきわめて大である。
（変形例）
本実施形態は前記した実施形態に限定されるものではなく、例えば以下のような種々の変形が可能である。
（１）前記実施形態は、複素展開係数のうち最大の絶対値を持つ係数の次数ｍを求め、次数ｍを境界とした展開係数間の自乗和の比を用いてドップラ速度を算出したが、複素展開係数間の自乗和の比を用いることもできる。
（２）前記実施形態は、超音波を用いて生体内部の血流などの速度分布や速度の空間分布を表示する医療診断用の血流計や血流描画装置に用いたが、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置に用いることができる。
本実施形態は前記した実施形態に限定されるものではなく、例えば以下のような種々の変形が可能である。
（１）前記実施形態は、複素展開係数のうち最大の絶対値を持つ係数の次数ｍを求め、次数ｍを境界とした展開係数間の自乗和の比を用いてドップラ速度を算出したが、複素展開係数間の自乗和の比を用いることもできる。
（２）前記実施形態は、超音波を用いて生体内部の血流などの速度分布や速度の空間分布を表示する医療診断用の血流計や血流描画装置に用いたが、電磁波を送受信することにより、雨雲など運動する反射物を検出して描画する気象レーダーや、飛行物体を検出する航空レーダー、あるいは、近づく物体を検出する衝突防止レーダーなどのパルス・ドップラ・レーダー装置に用いることができる。
Claims (9)
- 速度計測対象物に対してＮ個のパルス波を所定間隔で送信し、前記パルス波の受信エコー信号を送信時刻の順番に並べた時系列信号を用いて前記速度計測対象物を構成する運動反射物の速度を計測する速度計測方法であって、
前記Ｎ個の時系列信号を０次から（Ｎ−１）次までの離散的ルジャンドル関数を基底として展開する展開ステップと、
自然数ｎとするとき、前記展開ステップで展開した１次以上（Ｎ−１）次以下の（２ｎ−１）次の展開係数と（２ｎ＋１）次の展開係数との線形結合に虚数単位を乗じて、２ｎ次の展開係数と線形結合した２ｎ次の複素展開係数と、（２ｎ＋１）次の展開係数に虚数単位を乗じた後、２ｎ次の展開係数と（２ｎ＋２）次の展開係数と線形結合した（２ｎ＋１）次の複素展開係数とを算出する複素展開係数算出ステップと、
前記各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数ｍを決定する次数決定ステップと、
前記次数ｍに対応した前記各展開係数間あるいは前記各複素展開係数間の自乗和の比から前記運動反射物に関する符号付き速度信号を算出する速度信号算出ステップと
からなる速度計測方法。 - 前記展開ステップは、ｎ次ルジャンドル関数を行ベクトルＬ（ｎ）で表し、これをｎ＝０から（Ｎ−１）まで列方向に並べたＮ×Ｎ行列ＬＬを転置した転置行列tＬＬを生成する過程と、前記転置行列tＬＬに前記Ｎ×Ｎ行列ＬＬを前から乗算したものの逆行列（ＬＬ・tＬＬ)-1を算出する過程と、前記Ｎ×Ｎ行列ＬＬに前記逆行列（ＬＬ・tＬＬ)-1を前から乗算する過程とによって算出された行列ＦＦを要素が前記時系列信号である列ベクトルＳｄ（ｎ）に前から乗算するステップであることを特徴とする請求の範囲第１項に記載の速度計測方法。
- 静止反射物に対する前記受信エコーの大きさの変動あるいはドリフトの程度によって、前記複素展開係数を小さな次数から徐々に無視することを特徴とする請求の範囲第２項に記載の速度計測方法。
- 前記運動反射物の速度が有効速度検出範囲よりも小さい場合に検出速度ゼロを出力するか否かによって請求の範囲第３項に記載の速度計測方法の実行を判定するステップを備えたことを特徴とする速度計測方法。
- 速度計測対象物に対してＮ個のパルス波を所定間隔で送信し、前記パルス波の受信エコー信号を送信時刻の順番に並べた時系列信号を用いて前記速度計測対象物を構成する運動反射物の速度を計測する速度計測装置であって、
前記Ｎ個の時系列信号を０次から（Ｎ−１）次までの離散的ルジャンドル関数を基底として展開する展開手段と、
自然数ｎとするとき、前記展開ステップで展開した１次以上（Ｎ−１）次以下の（２ｎ−１）次の展開係数と（２ｎ＋１）次の展開係数との線形結合に虚数単位を乗じて、２ｎ次の展開係数と線形結合した２ｎ次の複素展開係数と、（２ｎ＋１）次の展開係数に虚数単位を乗じた後、２ｎ次の展開係数と（２ｎ＋２）次の展開係数と線形結合した（２ｎ＋１）次の複素展開係数とを算出する複素展開係数算出手段と、
前記各複素展開係数のうち最大の絶対値を持つ係数の自然数である次数ｍを決定する次数決定手段と、
前記次数ｍに対応した前記各展開係数間あるいは前記各複素展開係数間の自乗和の比から前記運動反射物に関する符号付き速度信号を算出する速度信号算出手段と
からなる速度計測装置。 - 前記展開手段は、ｎ次ルジャンドル関数を行ベクトルＬ（ｎ）で表し、これをｎ＝０から（Ｎ−１）まで列方向に並べたＮ×Ｎ行列ＬＬを転置した転置行列tＬＬを生成する過程と、前記転置行列tＬＬに前記Ｎ×Ｎ行列ＬＬを前から乗算したものの逆行列（ＬＬ・tＬＬ)-1を算出する過程と、前記Ｎ×Ｎ行列ＬＬに前記逆行列（ＬＬ・tＬＬ)-1を前から乗算する過程とによって算出された行列ＦＦを要素が前記時系列信号である列ベクトルＳｄ（ｎ）に前から乗算する手段であることを特徴とする請求の範囲第５項に記載の速度計測装置。
- 静止反射物に対する前記受信エコーの大きさの変動あるいはドリフトの程度によって、前記複素展開係数を小さな次数から徐々に無視することを特徴とする請求の範囲第６項に記載の速度計測装置。
- 前記パルス波は、超音波であり、
前記運動反射物は、血流であることを特徴とする請求の範囲第５項ないし第７項に記載の速度計測装置。 - 前記パルス波は、超音波であり、
前記速度計測対象物は、肝臓であり、
前記肝臓の速度の空間微分の分布及び血流像を表示する表示手段を備えたことを特徴とする請求の範囲第５項ないし７項に記載の速度計測装置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2007517722A JP5108512B2 (ja) | 2005-05-27 | 2006-01-06 | 速度計測方法及びこれを用いた速度計測装置 |
Applications Claiming Priority (4)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2005154844 | 2005-05-27 | ||
| JP2005154844 | 2005-05-27 | ||
| PCT/JP2006/300081 WO2006126303A1 (ja) | 2005-05-27 | 2006-01-06 | 速度計測方法及びこれを用いた速度計測装置 |
| JP2007517722A JP5108512B2 (ja) | 2005-05-27 | 2006-01-06 | 速度計測方法及びこれを用いた速度計測装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPWO2006126303A1 true JPWO2006126303A1 (ja) | 2008-12-25 |
| JP5108512B2 JP5108512B2 (ja) | 2012-12-26 |
Family
ID=37451736
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2007517722A Expired - Fee Related JP5108512B2 (ja) | 2005-05-27 | 2006-01-06 | 速度計測方法及びこれを用いた速度計測装置 |
Country Status (5)
| Country | Link |
|---|---|
| US (1) | US7946992B2 (ja) |
| EP (1) | EP1884196A4 (ja) |
| JP (1) | JP5108512B2 (ja) |
| CN (1) | CN101184443B (ja) |
| WO (1) | WO2006126303A1 (ja) |
Families Citing this family (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8585733B2 (en) | 2006-04-19 | 2013-11-19 | Vibrynt, Inc | Devices, tools and methods for performing minimally invasive abdominal surgical procedures |
| US9636076B2 (en) | 2011-03-29 | 2017-05-02 | Toshiba Medical Systems Corporation | X-ray CT apparatus and image processing method |
| JP5872323B2 (ja) * | 2011-03-29 | 2016-03-01 | 株式会社東芝 | Ｘ線ｃｔ装置及び画像処理方法 |
| JP6104920B2 (ja) | 2011-10-13 | 2017-03-29 | マシモ・コーポレイション | 医療用監視ハブ |
| CN102503138A (zh) * | 2011-10-18 | 2012-06-20 | 武汉理工大学 | 一种烧结法建筑装饰微晶玻璃的着色方法 |
| US20130172734A1 (en) * | 2011-12-30 | 2013-07-04 | General Electric Company | Flow measurement with time-resolved data |
| JP6253360B2 (ja) * | 2013-02-13 | 2017-12-27 | キヤノン株式会社 | 被検体情報取得装置、被検体情報取得方法、及びプログラム |
| TWI507709B (zh) * | 2013-03-20 | 2015-11-11 | Univ Nat Taiwan | 藉由格雷編碼激發之超音波都卜勒偵測方法 |
| KR20160046669A (ko) * | 2014-10-21 | 2016-04-29 | 알피니언메디칼시스템 주식회사 | 빔포밍 장치, 초음파 이미징 장치 및 빔포밍 방법 |
| KR102040853B1 (ko) * | 2015-03-27 | 2019-11-05 | 알피니언메디칼시스템 주식회사 | 공간 스무딩 연산이 간단한 빔포밍 장치, 초음파 이미징 장치 및 빔포밍 방법 |
| DE102015207894A1 (de) * | 2015-04-29 | 2016-11-03 | Bayer Pharma Aktiengesellschaft | Ermitteln der Geschwindigkeit eines Fluids mit Hilfe eines bildgebenden Verfahrens |
| CN112105952A (zh) * | 2018-05-09 | 2020-12-18 | 古野电气株式会社 | 气象雷达装置、气象观测方法及气象观测程序 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPS6480351A (en) * | 1987-09-24 | 1989-03-27 | Hitachi Medical Corp | Ultrasonic doppler meter |
| JPH0484953A (ja) * | 1990-07-30 | 1992-03-18 | Matsushita Electric Ind Co Ltd | 超音波ドプラ映像装置 |
| JP2003521341A (ja) * | 2000-01-31 | 2003-07-15 | アー．ヤー． アンゲルセン、ビョルン | 医療用超音波イメージングにおける位相面収差およびパルス残響の補正 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1003768B (zh) * | 1985-04-01 | 1989-04-05 | 复旦大学 | 一种双向定量测量血流绝对速度的超声多普勒方法与仪器 |
| US5228009A (en) | 1992-04-10 | 1993-07-13 | Diasonics, Inc. | Parametric clutter elimination |
-
2006
- 2006-01-06 EP EP06702059A patent/EP1884196A4/en not_active Withdrawn
- 2006-01-06 CN CN2006800186517A patent/CN101184443B/zh not_active Expired - Fee Related
- 2006-01-06 JP JP2007517722A patent/JP5108512B2/ja not_active Expired - Fee Related
- 2006-01-06 WO PCT/JP2006/300081 patent/WO2006126303A1/ja active Application Filing
- 2006-01-06 US US11/915,625 patent/US7946992B2/en not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPS6480351A (en) * | 1987-09-24 | 1989-03-27 | Hitachi Medical Corp | Ultrasonic doppler meter |
| JPH0484953A (ja) * | 1990-07-30 | 1992-03-18 | Matsushita Electric Ind Co Ltd | 超音波ドプラ映像装置 |
| JP2003521341A (ja) * | 2000-01-31 | 2003-07-15 | アー．ヤー． アンゲルセン、ビョルン | 医療用超音波イメージングにおける位相面収差およびパルス残響の補正 |
Non-Patent Citations (2)
| Title |
|---|
| CSNC200901611046, 複素ルジャンドル法によるクラッタに強いドプラ速度検出, "一般演題抄録（口演）その１", 超音波医学 Ｖｏｌ．３２， Ｓｕｐｐｌｅｍｅｎｔ ＪＡＰＡＮＥＳＥ ＪＯＵＲＮＡＬ ＯＦ ＭＥＤＩＣＡＬ ＵＬＴＲＡＳＯＮＩＣＳ, 2005, 第３２巻, ｐ．２８１, JP, （社）日本超音波医学会 * |
| JPN6009038695, Ｈａｎｓ Ｔｏｒｐ, "Ｃｌｕｔｔｅｒ Ｒｅｊｅｃｔｉｏｎ Ｆｉｌｔｅｒｓ ｉｎ Ｃｏｌｏｒ Ｆｌｏｗ Ｉｍａｇｉｎｇ： Ａ Ｔｈｅｏｒｅｔｉｃａｌ Ａｐｐｒｏａｃｈ", ＩＥＥＥ ＴＲＡＮＳＡＣＴＩＯＮＳ ＯＮ ＵＬＴＲＡＳＯＮＩＣＳ， ＦＥＲＲＯＥＬＥＣＴＲＩＣＳ， ＡＮＤ ＦＲＥＱＵＥＮＣＹ ＣＯＮＴＲＯＬ, 199703, ｖｏｌ．４４， ｎｏ．２, ｐｐ．４１７−４２４ * |
Also Published As
| Publication number | Publication date |
|---|---|
| EP1884196A4 (en) | 2009-12-23 |
| US20090177091A1 (en) | 2009-07-09 |
| CN101184443A (zh) | 2008-05-21 |
| JP5108512B2 (ja) | 2012-12-26 |
| EP1884196A1 (en) | 2008-02-06 |
| WO2006126303A1 (ja) | 2006-11-30 |
| CN101184443B (zh) | 2011-12-28 |
| US7946992B2 (en) | 2011-05-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP5108512B2 (ja) | 速度計測方法及びこれを用いた速度計測装置 | |
| Ebrahimkhanlou et al. | Damage localization in metallic plate structures using edge-reflected lamb waves | |
| Giusti et al. | ISAR image resolution enhancement: Compressive sensing versus state-of-the-art super-resolution techniques | |
| Aderson | Multi-dimensional velocity estimation with ultrasound using spatial quadrature | |
| Levine et al. | Model-based imaging of damage with Lamb waves via sparse reconstruction | |
| Stankovic et al. | Inverse radon transform–based micro-Doppler analysis from a reduced set of observations | |
| Zhao et al. | Mode identification and extraction of broadband ultrasonic guided waves | |
| JP4369427B2 (ja) | ドプラ速度検出装置及びそれを用いた超音波診断装置 | |
| Jia et al. | Rigid and elastic acoustic scattering signal separation for underwater target | |
| US11719613B2 (en) | Method and device for quantifying viscoelasticity of a medium | |
| CN104854473A (zh) | 使用多普勒聚焦的亚奈奎斯特雷达处理 | |
| Mozumi et al. | Anti-aliasing method for ultrasonic 2D phase-sensitive motion estimator | |
| Noumura et al. | Feasibility of low-frequency directive sound source with high range resolution using pulse compression technique | |
| CN107907591B (zh) | 多组分固液两相混合物组分浓度的超声检测系统和方法 | |
| Sando et al. | High-frequency ultrasonic airborne Doppler method for noncontact elasticity measurements of living tissues | |
| Tian et al. | Parametric modeling in food package defect imaging | |
| JPH05146438A (ja) | 超音波装置およびそれを用いたベクトルドプラ法 | |
| Xiaoming et al. | Application of Hilbert–Huang transform to laser Doppler velocimeter | |
| Liu et al. | Order analysis for measurement of Lamb wave propagation distance | |
| Ingram et al. | Ultrasonic array imaging through reverberating layers for industrial process analysis | |
| Latif et al. | Dispersion analysis of acoustic circumferential waves using time-frequency representations | |
| Nam et al. | Polar coordinate-based guided wave beamforming imaging using scanning LDV | |
| Kuga et al. | Alias-free estimation of blood velocity using phase difference and amplitude correlation of chirp echo spectrum | |
| Li | Enhanced Ultrasonic Imaging and Non-Destructive Evaluation Using Angle-of-Arrival Estimation in Dimension-Reduced Beam Space | |
| Tran et al. | A Detailed study of Sonar Tomographic Imaging. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20101109 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20110111 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20110308 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20110509 |
|
| A02 | Decision of refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A02 Effective date: 20120313 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20120611 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20120615 |
|
| A911 | Transfer to examiner for re-examination before appeal (zenchi) |
Free format text: JAPANESE INTERMEDIATE CODE: A911 Effective date: 20120705 |
|
| 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: 20121002 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20121005 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20151012 Year of fee payment: 3 |
|
| LAPS | Cancellation because of no payment of annual fees |