JP4566913B2 - 平均化反復法を用いたスパイクノイズ除去方法及びコンピュータプログラム - Google Patents
平均化反復法を用いたスパイクノイズ除去方法及びコンピュータプログラム Download PDFInfo
- Publication number
- JP4566913B2 JP4566913B2 JP2005515089A JP2005515089A JP4566913B2 JP 4566913 B2 JP4566913 B2 JP 4566913B2 JP 2005515089 A JP2005515089 A JP 2005515089A JP 2005515089 A JP2005515089 A JP 2005515089A JP 4566913 B2 JP4566913 B2 JP 4566913B2
- Authority
- JP
- Japan
- Prior art keywords
- value
- spike noise
- data
- predetermined
- point
- 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.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
- G01D3/00—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups
- G01D3/028—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups mitigating undesired influences, e.g. temperature, pressure
- G01D3/032—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups mitigating undesired influences, e.g. temperature, pressure affecting incoming signal, e.g. by averaging; gating undesired signals
Description
本発明は、データ波形に混入するスパイクノイズを、効果的に除去するスパイクノイズ除去方法及びスパイクノイズを効果的に除去するコンピュータプログラムに関する。
従来、物理測定、電磁測定、光学測定、ラマン分光測定等において、１次元、２次元又は３次元空間内でデータが計測されている。
計測により得られたデータ波形は、様々なノイズを含み、特に、近傍のデータ値に対して極端に大きな値あるいは小さな値のデータ（以下、スパイクノイズと称する）を含むことがある。スパイクノイズは、電気的衝撃、機械的衝撃、あるいはその他の外的要因により発生するが、データ波形の質を著しく低下させるため、除去すべき対象となる。スパイクノイズを除去する方法として、これまでに以下の単純平均化法、メディアンフィルター法、時系列法が提案されている。ここで、近傍とは、関心のあるデータ点から所定範囲内を言うものとする。
第一のスパイクノイズ除去方法として、特許第２６２９５１１号に記載された方法が挙げられる。この方法では、データ波形の１次微分値が閾値以上である場合にスパイクノイズを検出したと判断し、１次微分値の符号変化点をスパイクノイズ発生位置とみなす。次に、スパイクノイズ発生位置を挟む前後のデータ値の平均値で、スパイクノイズ発生位置のデータ値を１回だけ置換する。この方法を以下、単純平均化法と称する。
第二のスパイクノイズ除去方法として、メディアンフィルター法が挙げられる。この方法では、各データ点から所定範囲内にあるデータ点のデータ値の中央値で、前記各データ点の各データ値を置換する。この処理は、各データ点がスパイクノイズ点か否かとは無関係に行なわれ、処理対象にスパイクノイズ値を含むため、置換に用いる中央値が正常なデータ値であれば、スパイクノイズは効果的に除去される。この方法は２次元データ波形のほか、特に画像等の３次元データ波形において多用され、一般にメディアンフィルター法と呼ばれる。
第三のスパイクノイズ除去方法として、データの時間的相関性に着目する方法が挙げられる。高感度の電荷結合素子によりラマンスペクトルを計測するにあたり、データ波形に混入するスパイクノイズは、大気中を飛散する宇宙線が原因であることが知られており、その除去方法が、Ｈ．Ｔａｋｅｕｃｈｉ，Ｓ．Ｈａｓｈｉｍｏｔｏ，ａｎｄ Ｉ．Ｈａｒａｄａ，Ａｐｐｌ．Ｓｐｅｃｔｒｏｓｃ．，４７（１），１２９（１９９３）に記載されている。この方法では、データ波形を２回測定してデータ波形間の差を計算し、差が閾値以上となる箇所をスパイクノイズ発生位置とみなし、データ値の大きな方をスパイクノイズ、小さな方を真のデータと判断して、前者を後者で置換する。これが第三のスパイクノイズ除去方法であり、以下時系列法と称する。
計測により得られたデータ波形は、様々なノイズを含み、特に、近傍のデータ値に対して極端に大きな値あるいは小さな値のデータ（以下、スパイクノイズと称する）を含むことがある。スパイクノイズは、電気的衝撃、機械的衝撃、あるいはその他の外的要因により発生するが、データ波形の質を著しく低下させるため、除去すべき対象となる。スパイクノイズを除去する方法として、これまでに以下の単純平均化法、メディアンフィルター法、時系列法が提案されている。ここで、近傍とは、関心のあるデータ点から所定範囲内を言うものとする。
第一のスパイクノイズ除去方法として、特許第２６２９５１１号に記載された方法が挙げられる。この方法では、データ波形の１次微分値が閾値以上である場合にスパイクノイズを検出したと判断し、１次微分値の符号変化点をスパイクノイズ発生位置とみなす。次に、スパイクノイズ発生位置を挟む前後のデータ値の平均値で、スパイクノイズ発生位置のデータ値を１回だけ置換する。この方法を以下、単純平均化法と称する。
第二のスパイクノイズ除去方法として、メディアンフィルター法が挙げられる。この方法では、各データ点から所定範囲内にあるデータ点のデータ値の中央値で、前記各データ点の各データ値を置換する。この処理は、各データ点がスパイクノイズ点か否かとは無関係に行なわれ、処理対象にスパイクノイズ値を含むため、置換に用いる中央値が正常なデータ値であれば、スパイクノイズは効果的に除去される。この方法は２次元データ波形のほか、特に画像等の３次元データ波形において多用され、一般にメディアンフィルター法と呼ばれる。
第三のスパイクノイズ除去方法として、データの時間的相関性に着目する方法が挙げられる。高感度の電荷結合素子によりラマンスペクトルを計測するにあたり、データ波形に混入するスパイクノイズは、大気中を飛散する宇宙線が原因であることが知られており、その除去方法が、Ｈ．Ｔａｋｅｕｃｈｉ，Ｓ．Ｈａｓｈｉｍｏｔｏ，ａｎｄ Ｉ．Ｈａｒａｄａ，Ａｐｐｌ．Ｓｐｅｃｔｒｏｓｃ．，４７（１），１２９（１９９３）に記載されている。この方法では、データ波形を２回測定してデータ波形間の差を計算し、差が閾値以上となる箇所をスパイクノイズ発生位置とみなし、データ値の大きな方をスパイクノイズ、小さな方を真のデータと判断して、前者を後者で置換する。これが第三のスパイクノイズ除去方法であり、以下時系列法と称する。
前記の単純平均化法、メディアンフィルター法、時系列法は、いずれもスパイクノイズが少数で、 孤立して存在する場合は有効である。しかし実際には、スパイクノイズが多数混入したり２個以上連続することも多い。この場合、前記三つの方法では、実施例記載の通り、スパイクノイズが残存したり、元のデータ波形に歪みが発生する。スパイクノイズが多数混入したり、２個以上連続する場合に、生データ波形に歪みを発生させずにスパイクノイズを完全に除去することは、従来の技術では解決が困難であった。
本発明は、前述の実情に鑑みて提案されるものであって、２個以上連続したスパイクノイズを、データ波形に歪みを発生させずに、検出し除去するスパイクノイズ除去方法及びコンピュータプログラムを提供することを目的とする。
前記の課題は、スパイクノイズ範囲内の最大値あるいは最小値近傍の多数の点を除去対象とすること、除去対象となった個々の点の値を、除去対象データ近傍のデータの平均値で繰り返し置換することにより解決できる。このような、課題の解決手段となる本発明の方法を、平均化反復法と称する。
平均化反復法において、平均値の置換を繰り返すにあたり、除去対象のデータ点よりチャンネル番号の大きなデータ点のデータ値と、チャンネル番号の小さなデータ点のデータ値とを用いて平均値を計算する。これには、除去対象のデータ点よりチャンネル番号が１だけ大きいデータ点のデータ値と、チャンネル番号が１だけ小さいデータ点のデータ値と（２つの最隣接点）を用いて平均値を計算する場合も含まれる。
本発明に係るスパイクノイズ除去方法は、物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが正値でΔｙｇ，ｇ＋１が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最大値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）前記最大値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差が、当該所定閾値以上であるとき、最大値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、を有する。
または、本発明に係るスパイクノイズ除去方法は、物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが負値で｜Δｙｇ，ｇ＋１｜が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最小値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）前記最小値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、最小値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、を有する。
ステップ（３）で前記差の絶対値が、当該所定閾値以上であるとき、前記ステップ（１）（２）（３）の後に、（４）データ値ｙｍ−ｃ，ｙｍ＋ｄ（ただしｃ，ｄは１≦ａ＜ｃ≦Ｘ２，１≦ｂ＜ｄ≦Ｘ２を満たす特定の値、Ｘ２は所定の第２間隔）の平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２を計算するステップと、（５）前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２と前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２との差を計算するステップと、（６）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２で置換するステップと、（７）値ａを値ｃで置換し、値ｂを値ｄで置換するステップと、（８）ａ＜ｃ≦Ｘ２，ｂ＜ｄ≦Ｘ２の範囲の各ｃ，各ｄに対して、前記（４）（５）（６）（７）を行うステップ（Ｓ３３）と、をさらに有することが好ましい。
前記ステップ（１）（２）（３）（４）（５）（６）（７）（８）の後に、（９）前記スパイクノイズ点（ｍ，ｙｍ）から所定第３間隔Ｘ１−１の範囲内に位置する各データ点に対して、前記スパイクノイズ点（ｍ，ｙｍ）に対する、前記（１）（２）（３）（４）（５）（６）（７）（８）の処理と同様の処理を行い、当該スパイクノイズ範囲内のスパイクノイズが除去されたと判断するステップ（Ｓ３２）と、をさらに有することが好ましい。
Ｘ２≧Ｘ１−１であることが好ましい。ここで、Ｘ２はデータ列の最大数である。これは、ＣＣＤのような検出器のチャンネル数に相当し、例えば５１２，１０２４のような値である。
各測定点の位置は、受光素子のチャンネル位置であることが好ましい。
前記スパイクノイズ範囲の終点を判断するステップは、前記始点から所定の第１間隔Ｘ１の範囲内において、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、前記閾値（Ｙ１）以上にならないとき、当該範囲内に、スパイクノイズは存在しないとみなすステップ（Ｓ２５ａ）を有することが好ましい。
または、前記スパイクノイズ範囲の終点を判断するステップは、前記始点から所定の第１間隔Ｘ１の範囲内において、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上にならないとき、当該範囲内に、スパイクノイズは存在しないとみなすステップ（Ｓ２５ａ）を有することが好ましい。
本発明に係るコンピュータプログラムは、前述のスパイクノイズ除去方法の計算方法をコンピュータに実行させるものである。
本発明に係るコンピュータプログラムは、物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが正値でΔｙｇ，ｇ＋１が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最大値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）前記最大値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差が、当該所定閾値以上であるとき、最大値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、ステップ（３）で前記差が、当該所定閾値以上であるとき、前記ステップ（１）（２）（３）の後に、（４）データ値ｙｍ−ｃ，ｙｍ＋ｄ（ただしｃ，ｄは１≦ａ＜ｃ≦Ｘ２，１≦ｂ＜ｄ≦Ｘ２を満たす特定の値、Ｘ２は所定の第２間隔）の平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２を計算するステップと、（５）前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２と前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２との差を計算するステップと、（６）前記差が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差が、当該所定閾値以上であるとき、前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２で置換するステップと、（７）値ａを値ｃで置換し、値ｂを値ｄで置換するステップと、（８）ａ＜ｃ≦Ｘ２，ｂ＜ｄ≦Ｘ２の範囲の各ｃ，各ｄに対して、前記（４）（５）（６）（７）を行うステップ（Ｓ３３）と、（９）前記スパイクノイズ点（ｍ，ｙｍ）から所定第３間隔Ｘ１−１の範囲内に位置する各データ点に対して、前記スパイクノイズ点（ｍ，ｙｍ）に対する、前記（１）（２）（３）（４）（５）（６）（７）（８）の処理と同様の処理を行い、当該スパイクノイズ範囲内のスパイクノイズが除去されたと判断するステップ（Ｓ３２）と、をコンピュータに実行させる。
または、本発明に係るコンピュータプログラムは、物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが負値で｜Δｙｇ，ｇ＋１｜が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最小値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）前記最小値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、最小値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、ステップ（３）で前記差の絶対値が、当該所定閾値以上であるとき、前記ステップ（１）（２）（３）の後に、（４）データ値ｙｍ−ｃ，ｙｍ＋ｄ（ただしｃ，ｄは１≦ａ＜ｃ≦Ｘ２，１≦ｂ＜ｄ≦Ｘ２を満たす特定の値、Ｘ２は所定の第２間隔）の平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２を計算するステップと、（５）前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２と前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２との差を計算するステップと、（６）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２で置換するステップと、（７）値ａを値ｃで置換し、値ｂを値ｄで置換するステップと、（８）ａ＜ｃ≦Ｘ２，ｂ＜ｄ≦Ｘ２の範囲の各ｃ，各ｄに対して、前記（４）（５）（６）（７）を行うステップ（Ｓ３３）と、（９）前記スパイクノイズ点（ｍ，ｙｍ）から所定第３間隔Ｘ１−１の範囲内に位置する各データ点に対して、前記スパイクノイズ点（ｍ，ｙｍ）に対する、前記（１）（２）（３）（４）（５）（６）（７）（８）の処理と同様の処理を行い、当該スパイクノイズ範囲内のスパイクノイズが除去されたと判断するステップ（Ｓ３２）と、をコンピュータに実行させる。
本発明は、物理的測定、電磁測定、光学測定、ラマン分光測定に適用することが好ましい。また、本発明は、２次元又は３次元のデータの計測に適用することが好ましい。ここで、物理的測定とは、計測機器及び分析機器を用いて取得した信号等をいう。
本発明の方法（平均化反復法）によれば、データ波形中の２個以上連続したスパイクノイズを、データ波形に歪みを発生させずに、検出し除去することが可能である。
本発明は、前述の実情に鑑みて提案されるものであって、２個以上連続したスパイクノイズを、データ波形に歪みを発生させずに、検出し除去するスパイクノイズ除去方法及びコンピュータプログラムを提供することを目的とする。
前記の課題は、スパイクノイズ範囲内の最大値あるいは最小値近傍の多数の点を除去対象とすること、除去対象となった個々の点の値を、除去対象データ近傍のデータの平均値で繰り返し置換することにより解決できる。このような、課題の解決手段となる本発明の方法を、平均化反復法と称する。
平均化反復法において、平均値の置換を繰り返すにあたり、除去対象のデータ点よりチャンネル番号の大きなデータ点のデータ値と、チャンネル番号の小さなデータ点のデータ値とを用いて平均値を計算する。これには、除去対象のデータ点よりチャンネル番号が１だけ大きいデータ点のデータ値と、チャンネル番号が１だけ小さいデータ点のデータ値と（２つの最隣接点）を用いて平均値を計算する場合も含まれる。
本発明に係るスパイクノイズ除去方法は、物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが正値でΔｙｇ，ｇ＋１が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最大値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）前記最大値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差が、当該所定閾値以上であるとき、最大値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、を有する。
または、本発明に係るスパイクノイズ除去方法は、物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが負値で｜Δｙｇ，ｇ＋１｜が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最小値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）前記最小値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、最小値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、を有する。
ステップ（３）で前記差の絶対値が、当該所定閾値以上であるとき、前記ステップ（１）（２）（３）の後に、（４）データ値ｙｍ−ｃ，ｙｍ＋ｄ（ただしｃ，ｄは１≦ａ＜ｃ≦Ｘ２，１≦ｂ＜ｄ≦Ｘ２を満たす特定の値、Ｘ２は所定の第２間隔）の平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２を計算するステップと、（５）前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２と前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２との差を計算するステップと、（６）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２で置換するステップと、（７）値ａを値ｃで置換し、値ｂを値ｄで置換するステップと、（８）ａ＜ｃ≦Ｘ２，ｂ＜ｄ≦Ｘ２の範囲の各ｃ，各ｄに対して、前記（４）（５）（６）（７）を行うステップ（Ｓ３３）と、をさらに有することが好ましい。
前記ステップ（１）（２）（３）（４）（５）（６）（７）（８）の後に、（９）前記スパイクノイズ点（ｍ，ｙｍ）から所定第３間隔Ｘ１−１の範囲内に位置する各データ点に対して、前記スパイクノイズ点（ｍ，ｙｍ）に対する、前記（１）（２）（３）（４）（５）（６）（７）（８）の処理と同様の処理を行い、当該スパイクノイズ範囲内のスパイクノイズが除去されたと判断するステップ（Ｓ３２）と、をさらに有することが好ましい。
Ｘ２≧Ｘ１−１であることが好ましい。ここで、Ｘ２はデータ列の最大数である。これは、ＣＣＤのような検出器のチャンネル数に相当し、例えば５１２，１０２４のような値である。
各測定点の位置は、受光素子のチャンネル位置であることが好ましい。
前記スパイクノイズ範囲の終点を判断するステップは、前記始点から所定の第１間隔Ｘ１の範囲内において、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、前記閾値（Ｙ１）以上にならないとき、当該範囲内に、スパイクノイズは存在しないとみなすステップ（Ｓ２５ａ）を有することが好ましい。
または、前記スパイクノイズ範囲の終点を判断するステップは、前記始点から所定の第１間隔Ｘ１の範囲内において、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上にならないとき、当該範囲内に、スパイクノイズは存在しないとみなすステップ（Ｓ２５ａ）を有することが好ましい。
本発明に係るコンピュータプログラムは、前述のスパイクノイズ除去方法の計算方法をコンピュータに実行させるものである。
本発明に係るコンピュータプログラムは、物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが正値でΔｙｇ，ｇ＋１が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最大値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）前記最大値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差が、当該所定閾値以上であるとき、最大値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、ステップ（３）で前記差が、当該所定閾値以上であるとき、前記ステップ（１）（２）（３）の後に、（４）データ値ｙｍ−ｃ，ｙｍ＋ｄ（ただしｃ，ｄは１≦ａ＜ｃ≦Ｘ２，１≦ｂ＜ｄ≦Ｘ２を満たす特定の値、Ｘ２は所定の第２間隔）の平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２を計算するステップと、（５）前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２と前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２との差を計算するステップと、（６）前記差が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差が、当該所定閾値以上であるとき、前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２で置換するステップと、（７）値ａを値ｃで置換し、値ｂを値ｄで置換するステップと、（８）ａ＜ｃ≦Ｘ２，ｂ＜ｄ≦Ｘ２の範囲の各ｃ，各ｄに対して、前記（４）（５）（６）（７）を行うステップ（Ｓ３３）と、（９）前記スパイクノイズ点（ｍ，ｙｍ）から所定第３間隔Ｘ１−１の範囲内に位置する各データ点に対して、前記スパイクノイズ点（ｍ，ｙｍ）に対する、前記（１）（２）（３）（４）（５）（６）（７）（８）の処理と同様の処理を行い、当該スパイクノイズ範囲内のスパイクノイズが除去されたと判断するステップ（Ｓ３２）と、をコンピュータに実行させる。
または、本発明に係るコンピュータプログラムは、物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが負値で｜Δｙｇ，ｇ＋１｜が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最小値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）前記最小値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、最小値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、ステップ（３）で前記差の絶対値が、当該所定閾値以上であるとき、前記ステップ（１）（２）（３）の後に、（４）データ値ｙｍ−ｃ，ｙｍ＋ｄ（ただしｃ，ｄは１≦ａ＜ｃ≦Ｘ２，１≦ｂ＜ｄ≦Ｘ２を満たす特定の値、Ｘ２は所定の第２間隔）の平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２を計算するステップと、（５）前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２と前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２との差を計算するステップと、（６）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２で置換するステップと、（７）値ａを値ｃで置換し、値ｂを値ｄで置換するステップと、（８）ａ＜ｃ≦Ｘ２，ｂ＜ｄ≦Ｘ２の範囲の各ｃ，各ｄに対して、前記（４）（５）（６）（７）を行うステップ（Ｓ３３）と、（９）前記スパイクノイズ点（ｍ，ｙｍ）から所定第３間隔Ｘ１−１の範囲内に位置する各データ点に対して、前記スパイクノイズ点（ｍ，ｙｍ）に対する、前記（１）（２）（３）（４）（５）（６）（７）（８）の処理と同様の処理を行い、当該スパイクノイズ範囲内のスパイクノイズが除去されたと判断するステップ（Ｓ３２）と、をコンピュータに実行させる。
本発明は、物理的測定、電磁測定、光学測定、ラマン分光測定に適用することが好ましい。また、本発明は、２次元又は３次元のデータの計測に適用することが好ましい。ここで、物理的測定とは、計測機器及び分析機器を用いて取得した信号等をいう。
本発明の方法（平均化反復法）によれば、データ波形中の２個以上連続したスパイクノイズを、データ波形に歪みを発生させずに、検出し除去することが可能である。
図１Ａ及び１Ｂは、平均値による置換の反復で、スパイクノイズを除去してゆくプロセスを、データ列を例に示した図である。
図２は、本実施形態の計算方法（平均化反復法）を説明するフローチャートである。
図３は、閾値の決定の過程を説明するフローチャートである。
図４は、スパイクノイズの検出とスパイクノイズ範囲の決定の過程を説明するフローチャートである。
図５は、スパイクノイズの除去の過程を説明するフローチャートである。
図６Ａ及び６Ｂは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
図６Ｃ及び６Ｄは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
図６Ｅ及び６Ｆは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
図６Ｇ及び６Ｈは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
図６Ｉは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
図２は、本実施形態の計算方法（平均化反復法）を説明するフローチャートである。
図３は、閾値の決定の過程を説明するフローチャートである。
図４は、スパイクノイズの検出とスパイクノイズ範囲の決定の過程を説明するフローチャートである。
図５は、スパイクノイズの除去の過程を説明するフローチャートである。
図６Ａ及び６Ｂは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
図６Ｃ及び６Ｄは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
図６Ｅ及び６Ｆは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
図６Ｇ及び６Ｈは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
図６Ｉは、アクリル板のラマンスペクトルと、各種スパイクノイズ除去方法の適用結果を示した図である。
以下、図面を参照して、本発明の実施形態を説明する。
本実施形態は、スパイクノイズ検出・除去に必要な閾値の決定、スパイクノイズの検出とスパイクノイズ範囲の決定、スパイクノイズの除去という３つの過程（ステップ）からなる。
より詳細には、本実施形態のスパイクノイズ除去方法は、物理的測定により得られた各点の上のデータについて、データ波形（ｉ，ｙｉ）の（ａ）１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが正値でΔｙｇ，ｇ＋１が、又は（ｂ）Δｙｇ，ｇ＋１が負値で｜Δｙｇ，ｇ＋１｜が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、（ａ）１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、又は（ｂ）Δｙｈ，ｈ＋１が正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における（ａ）最大値（ｙｍ）又は（ｂ）最小値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）（ａ）最大値（ｙｍ）又は（ｂ）最小値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、（ａ）最大値（ｙｍ）又は（ｂ）最小値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、（４）データ値ｙｍ−ｃ，ｙｍ＋ｄ（ただしｃ，ｄは１≦ａ＜ｃ≦Ｘ２，１≦ｂ＜ｄ≦Ｘ２を満たす特定の値、Ｘ２は所定の第２間隔）の平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２を計算するステップと、（５）前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２と前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２との差を計算するステップと、（６）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２で置換するステップと、（７）文字ａの値を文字ｃの値で置換し、文字ｂの値を文字ｄの値で置換するステップと、（８）ａ＜ｃ≦Ｘ２，ｂ＜ｄ≦Ｘ２の範囲の各ｃ，各ｄに対して、前記（４）（５）（６）（７）を行い、前記スパイクノイズ点（ｍ，ｙｍ）に対する処理が終了したと判断するステップ（Ｓ３３）と、（９）前記スパイクノイズ点（ｍ，ｙｍ）から所定第３間隔Ｘ１−１の範囲内に位置する各データ点に対して、前記スパイクノイズ点（ｍ，ｙｍ）に対する、前記（１）（２）（３）（４）（５）（６）（７）（８）の処理と同様の処理を行い、当該スパイクノイズ範囲内のスパイクノイズが除去されたと判断するステップ（Ｓ３２）とを有する。
図１Ａ及び１Ｂは、スパイクノイズの除去の過程の実施形態の考え方を示す。
図１Ａは、スパイクノイズの除去の操作の説明のために作成した２次元データ波形であり、スパイクノイズ範囲を４≦ｉ≦７とし、ｙ５＝５０，ｙ６＝４６を、バックグラウンドに混入した２本の連続したスパイクノイズと仮定する。図中折線（ａ１）は生データ、折線（ａ２）はスパイクノイズ除去後のデータを示す。
また、前記閾値Ｘ１−１を３、計算終了判定のための閾値（前記閾値Ｙ２）として、ＮＬ（ノイズレベル）＝４と仮定し、図１Ｂには、スパイクノイズを含む生データが、スパイクノイズの除去過程で変化してゆく様子を示した。
図１Ｂ第２列において、生データのスパイクノイズ範囲４≦ｉ≦７内の最大値はｙ５＝５０であり、これが最初の除去対象となりうる。ｙ５＝５０と平均値ｙ５１＝（ｙ４＋ｙ６）／２＝２５との差を計算すると｜ｙ５−ｙ５１｜＝２５となり、この値は（ノイズレベル）＝４よりも大きい。そこで、ｙ５＝５０をｙ５１＝２５で置換する。次に、ｙ５１＝２５と平均値ｙ５２＝（ｙ３＋ｙ７）／２＝２との差を計算すると｜ｙ５２−ｙ５１｜＝２３となり、この値は（ノイズレベル）＝４よりも大きい。そこで、ｙ５１＝２５をｙ５２＝２で置換する。さらに、ｙ５２＝２と平均値ｙ５３＝（ｙ２＋ｙ８）／２＝３との差を計算すると｜ｙ５３−ｙ５２｜＝１となり、この値は（ノイズレベル）＝４よりも小さい。そこで、スパイクノイズｙ５＝５０の除去が終了したと判断し、ｙ５２＝２で置換した状態で計算を終了する。この段階でのデータ波形は、図１Ｂ第３列に示した通りである。上記と同様の処理を、ｉ＝５，５−１，５＋１，５−２，５＋２，５−３，５＋３＝５，４，６，３，７，２，８なるｙｉについて行う（図５のステップＳ２２，Ｓ３２において、Ｘ１＝４とした場合に相当する）。ｙ５に関する処理は終了したので、ｙ５に隣接するｙ４に関する処理を説明する。
ｙ４＝５とｙ４１＝（ｙ３＋ｙ５）／２＝２．５との差を計算すると、｜ｙ４−ｙ４１｜＝２．５となり、この値は（ノイズレベル）＝４よりも小さい。そこで、ｙ４＝５をスパイクノイズとはみなさずに、ｙ４＝５のまま計算を終了する。これにより、データ波形は図１Ｂ第３列の状態を保持し、図１Ｂ第４列のようになる。次に、ｙ４と同様にｙ５に隣接するｙ６に関する処理を説明する。
ｙ６＝４５と平均値ｙ６１＝（ｙ５＋ｙ７）／２＝１．５との差を計算すると｜ｙ６−ｙ６１｜＝４３．５となり、この値は（ノイズレベル）＝４よりも大きい。そこで、ｙ６＝４５をｙ６１＝１．５で置換する。次に、ｙ６１＝１．５と平均値ｙ６２＝（ｙ４＋ｙ８）／２＝５との差を計算すると｜ｙ６２−ｙ６１｜＝３．５となり、この値は（ノイズレベル）＝４よりも小さい。そこで、スパイクノイズｙ６＝４５の除去が終了したと判断し、ｙ６１＝１．５で置換した状態で計算を終了する。これにより、データ波形は、図１Ｂ第５列のように変化する。
上記ｙ５，ｙ４，ｙ６に関する処理と同様の処理を、ｙ３，ｙ７，ｙ２，ｙ８にも施すと、ｙ４の場合と同様に個々の値は変化せず、最終的なデータ波形は図１Ｂ第６列のようになる。以上の過程により、生データにおいて連続する２本のスパイクノイズｙ５＝５０，ｙ６＝４５が完全に除去され、平均値による置換を、範囲を拡大しながら反復することの利点が示された。
なお、単純平均化法（特許第２６２９５１１号の方法）では、生データは図１Ｂ第７列のように変化する。特許第２６２９５１１号の方法では、スパイクノイズ範囲は５≦ｉ≦６となり、スパイクノイズ範囲を挟むｙ４＝５，ｙ７＝１の平均値３でスパイクノイズ範囲内の値が全て置換されるためである。このように単純平均化法では、一部のデータを全て同一の値で置換するため、波形に歪みを発生することが予想される。また、メディアンフィルター法で３チャンネル内の中央値による置換を仮定すると、生データのｙ０，ｙ１，ｙ２，ｙ３，ｙ４，ｙ５，ｙ６，ｙ７，ｙ８は順に０，１，３，３，５，４５，４５，５，４，４，２に変化し、スパイクノイズが残存することになる。単純平均化法で発生が予想される波形の歪みは、図１のようにバックグラウンド上ではなく、シグナル波形上にスパイクノイズが混入した場合に顕著である。具体例として、特許第２６２９５１１号図８の測定時間３．２秒付近の波形に観測されるショルダーピークが挙げられる。本実施形態では、スパイクノイズ近傍のシグナルデータ値の平均値で置換を行う。このため置換後の波形には、スパイクノイズ近傍のシグナルデータの特徴が反映され、単純平均化法、メディアンフィルター法の場合よりも、置換後の波形の歪みが小さくなる。
図２〜図５は、本発明の実施形態を示すフローチャートである。図２に示すように、ステップＳ１０でスパイクノイズの検出・除去で必要となる閾値の決定を行い、ステップＳ２０でスパイクノイズの検出とスパイクノイズ範囲の決定を行い、ステップＳ３０はスパイクノイズの除去を行う。
後述するように、オペレータは、パラメータとして、ステップＳ１０で、バックグラウンド上の２点の位置、および認識パラメータ（ＲｅｃｏｇＰＡ）、およびノイズ比パラメータ（ＲａｔｉｏＰＡ）を決定し、Ｓ２０でスパイクノイズレンジ（ＣＲＲ）を決定する。
バックグラウンド上の上記の２点の回帰直線とデータの差の標準偏差がＳＩＧとして算出され、前記閾値Ｙ１に相当するのはＮＴ＝ＳＩＧ×ＲｅｃｏｇＰＡとなり、前記閾値Ｘ１−１に相当するのはＣＲＲとなり、前記閾値Ｙ２に相当するのはＮＬ＝ＳＩＧ×ＲａｔｉｏＰＡとなる。なお、本発明の実施形態はこのような閾値の決定法に限定されるものではなく、別の方法により前記閾値Ｙ１およびＹ２を指定することも可能である。
ステップＳ１０、Ｓ２０でパラメータが、一度指定されると、該パラメータはステップＳ２０、Ｓ３０で利用され、計算は自動的に進行する。
なお、図２〜図５において、オペレータによる指定が含まれる過程は平行四辺形で表示され、それ以外の過程は長方形あるいはひし形で表示されている。なお、ひし形は条件による分岐を意味している。また、連続するデータ間での１次微分値を差分と定義し、差分を用いて計算過程を説明してある。
以下、図３〜図５のフローチャートを図１のスパイクノイズ除去過程を参照して説明する。図３のステップＳ１１〜Ｓ１６は図２のステップＳ１０を詳細に示し、図４のステップＳ２１〜Ｓ２６、Ｓ２５ａは図２のステップＳ２０を詳細に示し、図５のステップＳ３１〜Ｓ３５、Ｓ３５ａは図２のステップＳ３０を詳細に示す。
ここで、前記閾値をＹ１＝１２（ＳＩＧ×ＲｅｃｏｇＰＡ＝１２），Ｘ１＝４（ＣＲＲ＝３），Ｙ２＝４（ＳＩＧ×ＲａｔｉｏＰＡ＝４、図１Ａにおけるノイズレベル），Ｘ２＝６（２×ＣＲＲ＝６）と仮定する。後述するように、これらの値は、図３のステップＳ１５、Ｓ１６、図４のステップＳ２２で指定される。
まず、図３に示すスパイクノイズの検出・除去で必要となる閾値の決定について説明する。
ステップＳ１１では、バックグラウンド上の２点を指定し、ステップＳ１２では、２点間の回帰直線とデータ値の差の標準偏差を計算してＳＩＧとする。ステップＳ１３で認識パラメータＲｅｃｏｇＰＡを指定し、ステップＳ１４でノイズ比パラメータＲａｔｉｏＰＡを指定する。ステップＳ１５では、ＮＴ＝ＳＩＧ×ＲｅｃｏｇＰＡを閾値（Ｙ１）として算出し、ステップＳ１６では、ＮＴ＝ＳＩＧ×ＲａｔｉｏＰＡを閾値（Ｙ２）として算出する。
次に、図４に示すスパイクノイズの検出とスパイクノイズ範囲の決定について説明する。
ステップＳ２１では、検出するスパイクノイズピーク値の正負を指定する（以下、スパイクノイズピーク値が正の場合を仮定する）。ステップＳ２２では、検出・除去の対象として許容するスパイクノイズ幅のチャンネル数の最大値より、２減じた値を、スパイクノイズレンジ（ＣＲＲ）として指定する。ここで、ＣＲＲ＝Ｘ１−１とする。ステップＳ２３では、データ（ｉ，ｙｉ）の差分Δｙｉ，ｉ＋１＝ｙｉ＋１−ｙｉを計算する。
パラメータの指定および計算がステップＳ２３まで進行した後、ステップＳ２４において、Δｙｊ，ｊ＋１＝ｙｊ＋１−ｙｊ＞＋ＮＴのとき、ｊをスパイクノイズの始まるチャンネルとする。図１Ａ及び１Ｂの例の場合、Δｙ４，５＝ｙ５−ｙ４＝４５＞１２であるのでｊ＝４である。
ステップＳ２５においてΔｙｊ，ｊ＋ｑ＋１＝ｙｊ＋ｑ＋１−ｙｊ＜−ＮＴのとき、ｊ＋ｑ＋１をスパイクノイズが終了するチャンネルとする。ただし、１≦ｑ≦ＣＲＲとする。図１の例の場合、Δｙ６，７＝ｙ７−ｙ６＝−４４＜−１２であり、ｊ＋ｑ＋１＝７であるから、ｑ＝２となる。
ステップＳ２６において、スパイクノイズ範囲はｊ≦ｉ≦ｊ＋ｑ＋１であるから、スパイクノイズ幅はチャンネル数ｑ＋２（間隔ｑ＋１）である。ただし、１≦ｑ≦ＣＲＲである。図１の例の場合、スパイクノイズ範囲は４≦ｉ≦７、ｑ＝２となる。ここでＣＲＲ＝３であるから、１≦ｑ≦ＣＲＲが成立する。ここで、スパイクノイズ範囲は４〜７チャンネルであり、スパイクノイズ幅は４チャンネルである。
ステップＳ２５ａでは、１≦ｑ≦ＣＲＲの範囲で差分がＳ２５の条件を満たさないときは、スパイクノイズは存在しないとみなし、計算を終了する。このステップＳ２５ａの過程は、スパイクノイズ幅が過剰に大きくなることを防止するために設けた。
ステップＳ２５ａの過程を省略すると、スパイクノイズ幅が、正常な信号ピークのバンド幅に相当する値を持つ場合、正常な信号ピークを広範囲におよぶスパイクノイズと判断し、このピークを図５のステップＳ３１以下の過程で除去してしまう可能性がある。
図１Ａ及び１Ｂの例の場合、ＣＲＲ＝３であるからスパイクノイズ幅のチャンネル数ｑ＋２は５以下となる必要がある。実際にはｑ＋２＝４であり、上記制限内であることがわかる。このようにＣＲＲの設定は、スパイクノイズ幅を制限し、連続したスパイクノイズと、正常な信号ピークとを区別する機能を持つ。
次に、図５に示すスパイクノイズの除去について説明する。
ステップＳ３１において、スパイクノイズ範囲ｊ≦ｉ≦ｊ＋ｑ＋１内でｙｉが最大の点（ｐ，ｙｐ）を特定する。図１の例の場合、スパイクノイズ範囲４≦ｉ≦７内でｙｉが最大の点は（５，ｙ５）であるので、ｐ＝５である。
ステップＳ３２では、ステップＳ３３〜Ｓ３５、Ｓ３５ａをｉ＝ｐ，ｐ−１，ｐ＋１，‥，ｐ−ＣＲＲ，ｐ＋ＣＲＲについて実行する。
ステップＳ３３では、ステップＳ３４、Ｓ３５、Ｓ３５ａを０≦ｒ≦２×ＣＲＲについて実行する。ここで、２×ＣＲＲ＝Ｘ２である。図１の例の場合、前記（５，ｙ５）を処理するにあたり、ｒは２×ＣＲＲ＝６以下となる必要がある。本例では、前述の通り、ｒ＝０，１，２で（５，ｙ５）の処理が完了しており、これは上記制限内であることがわかる。
ステップＳ３４では、ｙｉｒ＝（ｙｉ−ｒ＋ｙｉ＋ｒ）／２を計算する。
ステップＳ３４の後処理は分岐し、｜ｙｉｒ＋１−ｙｉｒ｜＞ＮＬの場合に処理はステップＳ３５に進み、｜ｙｉｒ＋１−ｙｉｒ｜＜ＮＬの場合に処理はステップＳ３５ａに進む。
このようにして（５，ｙ５）の処理がステップＳ３５ａで完了した後、再度ステップＳ３２に戻ることになる。図１Ａ及び１Ｂの例の場合、ステップＳ３２〜Ｓ３５ａで、上記と同様の処理を、ｉ＝５，５−１，５＋１，５−２，５＋２，５−３，５＋３＝５，４，６，３，７，２，８なるｙｉに施すことになる。そこで、ｙ５に関する説明の次に、ｙ５に隣接するｙ４に関する説明を行う。
ステップＳ３３〜Ｓ３５ａで、（４，ｙ４）を処理するにあたり、ｒは２×ＣＲＲ＝６以下となる必要がある。本例では、前述の通り、ｒ＝０で（４，ｙ４）の処理が完了しており、これは上記制限内であることがわかる。このようにして（４，ｙ４）の処理がステップＳ３５ａで完了し、再度ステップＳ３２に戻ることになる。次にｙ４と同様にｙ５に隣接するｙ６の処理に関する説明を行う。
ステップＳ３３〜Ｓ３５ａで、（６，ｙ６）を処理するにあたり、ｒは２×ＣＲＲ＝６以下となる必要がある。本例では、前述の通り、ｒ＝０，１で（６，ｙ６）の処理が完了しており、これは上記制限内であることがわかる。このようにして（６，ｙ６）の処理がステップＳ３５ａで完了した後、再度ステップＳ３２に戻ることになる。
ｙ５，ｙ４，ｙ６に関する処理と同様の処理をｙ３，ｙ７，ｙ２，ｙ８にも施すと、前述の通り、ｙ４の場合と同様に個々の値は変化せずｒ＝０で処理は完了し、ｒが２×ＣＲＲ＝６以下という制限内に収まっていることがわかる。このようにして、図１に示したように、連続した２本のスパイクノイズが除去されることになる。
本実施形態の本質は、平均値による置換の反復にあるが、これには二つの意味がある。第一は、図５のステップＳ３２の過程で、スパイクノイズ範囲内の最大値のほか、最大値近傍の複数の点が除去処理対象となること、第二は、図５のステップＳ３３の過程で、除去処理対象となった個々の点の値が、当該点周辺のベースラインに近似できるまで置換を繰り返すことである。
前記閾値Ｙ１，Ｘ１−１，Ｙ２，Ｘ２は、これまでの記述通り、ＳＩＧ×ＲｅｃｏｇＰＡ，ＣＲＲ，ＳＩＧ×ＲａｔｉｏＰＡ，２×ＣＲＲとするほか、別の方法でも指定できる。一例として、Ｘ２を１．５×ＣＲＲ，３×ＣＲＲ，４×ＣＲＲなどに設定してもよい。前記閾値Ｙ１，Ｘ１−１，Ｙ２，Ｘ２の指定法、あるいは上記のＸ２の指定法では、Ｘ２≧Ｘ１−１となる。特定のデータ値を平均値の置換で小さくする（ステップＳ３３）にあたり、Ｘ２≧Ｘ１−１は、重要な意味を持つ。Ｘ２≧Ｘ１−１の持つ意味は、上記の特定のデータ値から許容されるスパイクノイズの幅のチャンネル数（より２を減じた値）以上の距離内のデータを用いて、平均値を算出することにある。許容されるスパイクノイズの幅のチャンネル数（より２を減じた値）を下回る距離内のデータのみを用い、平均値を算出しても、スパイクノイズ範囲内のデータ値の平均値を利用することになり、スパイクノイズの除去が不十分になる可能性が高い。従って、前述の方法とは異なる方法で閾値Ｙ１，Ｘ１−１，Ｙ２，Ｘ２を指定する場合も、Ｘ２≧Ｘ１−１となることが望ましい。
閾値Ｘ１−１（ＣＲＲ）は、ステップＳ２２の通り、スパイクノイズ幅のチャンネル数の最大値より２減じた値を意味する。スパイクノイズ幅のチャンネル数の最大値がＸ１＋１となるから、スパイクノイズ範囲の始点から終点に至る間隔の最大値は、Ｘ１となる。Ｘ１は、請求項１記載の第１間隔に相当する。また、閾値Ｘ２（２×ＣＲＲ）は、ステップＳ３３の通り、除去処理対象となる特定のスパイクノイズのチャンネル番号と、平均値の算出に使用するデータのチャンネル番号との差の最大値を意味する。前記のデータのチャンネルと、前記のスパイクノイズのチャンネルとの差の最大値は、間隔にしてＸ２となる。Ｘ２は、請求項１、請求項２記載の第２間隔に相当する。さらにステップＳ３２の通り、除去処理対象となるデータ値のチャンネル番号と、スパイクノイズ範囲内の最大値のチャンネル番号との差の最大値は、Ｘ１−１（ＣＲＲ）である。これは、間隔にしてＸ１−１となり、請求項３記載の第３間隔に相当する。
前述の方法で閾値Ｙ１（ＳＩＧ×ＲｅｃｏｇＰＡ），Ｘ１−１（ＣＲＲ），Ｙ２（ＳＩＧ×ＲａｔｉｏＰＡ），Ｘ２（２×ＣＲＲ）を指定するにあたり、除去処理後も連続したスパイクノイズが残存するときはＸ１−１（ＣＲＲ）の値を大きくし、データ波形が歪むときはＸ１−１（ＣＲＲ）の値を小さくすればよい。閾値Ｘ１−１（ＣＲＲ）の値を大きくすること、ステップＳ３２が示すように、より多数のデータ値が除去処理の対象となる。また同時に閾値Ｘ２（２×ＣＲＲ）の値も大きくなり、ステップＳ３３が示すように、平均値の算出でより広範囲のデータ値を利用できるので、除去処理対象となった特定のデータ値を、スパイクノイズ範囲のデータ値の影響を受けずに、当該点周辺のベースラインに近似した値に置換できる。このように、閾値Ｘ１−１（ＣＲＲ）の値を大きくすることで、スパイクノイズのうち、高い連続性を示すもののを効果的に除去できる。上記の影響を考慮すれば、閾値Ｘ１−１（ＣＲＲ）は四つの閾値の中で最も本質的である。一方でステップＳ３２に代わり、ｉ＝ｐ，ｐ−１，ｐ＋１，…，ｐ−Ｘ３，ｐ＋Ｘ３とし、スパイクノイズの連続性が高いとき、Ｘ１−１，Ｘ２とは独立にＸ３を大きな値に設定して計算を実施してもよい。
また、滑らかなバックグラウンド上のスパイクノイズを除去するにあたり、除去処理後もスパイクノイズが残存するときは、ＲｅｃｏｇＰＡの値を大きくするのがよい。その理由は、滑らかなバックグラウンド上では標準偏差ＳＩＧの値が小さいために閾値Ｙ１（ＳＩＧ×ＲｅｃｏｇＰＡ）の値が小さくなり、ステップＳ２４、Ｓ２５でのスパイクノイズ範囲の決定が困難となるためである。
前述の方法で閾値Ｙ１，Ｘ１−１，Ｙ２，Ｘ２を指定するにあたり、連続したスパイクノイズが、滑らかなバックグラウンド上に存在する場合も含めて完全に除去され、かつデータ波形に歪みが発生しないよう各閾値を指定する必要があるが、前述の通り、前記の指定は過度の実験を必要としない。
以下、図６Ａ〜６Ｉに示した実施例により本実施形態を説明する。図６Ａは、ポリメタクリル酸メチルからなる板（以下、アクリル板と称する）のラマンスペクトルを、電荷結合素子を用い、スパイクノイズ除去機能なし、計測時間６００秒で計測した結果であり、宇宙線スパイクノイズが混入したスペクトルとなっている。図６Ｂ，６Ｄ，６Ｅ，６Ｆは、図６Ａに対して各種のスパイクノイズ除去方法を適用した結果である。図６Ｈは測定と同時に、時系列法によるスパイクノイズ除去機能を動作させ、スパイクノイズ除去を試みた結果である。図６Ｃ，６Ｇ，６Ｉは、それぞれ図６Ｂ，６Ｆ，６Ｈの一部拡大図である。ラマンスペクトルは、試料にレーザー光を照射し、試料が発する散乱光のうち照射光よりも低エネルギーのものについて、エネルギーの損失量毎に、散乱光の強度をプロットした２次元データである。エネルギーの損失量は、化学結合の種類や状態に依存するため、ここで得られたラマンスペクトルは、アクリル板固有のスペクトルとなる。
図６Ａ及び６Ｈのデータは、ダイオード励起可視レーザー（Ｓｐｅｃｔｒａ Ｐｈｙｓｉｃｓ Ｍｉｌｌｅｎｉａ）から発振する波長５３２ｎｍのレーザー光（強度０．６ｍＷ，ビーム径約０．２ｍｍ）を、厚さ２ｍｍのアクリル板に、アクリル板法線から約４５°の角度で照射して測定した。アクリル板から発生するラマン散乱光とレイリー散乱光を、アクリル板法線方向に設置したＦ／１の集光レンズで集光して平行光とした後、ノッチフィルター（Ｋａｉｓｅｒ Ｈｏｌｏｇｒａｐｈｉｃ Ｎｏｔｃｈ Ｆｉｌｔｅｒ ＨＮＰＦ−５３２．０−１．５）を通過させてレイリー散乱光をカットした。その後、Ｆ／７．５の集光レンズでシングルポリクロメータ（Ｊｏｂｉｎ Ｙｖｏｎ Ｒａｍａｎｏｒ Ｔ−６４０００，Ｆ／７．５）のスリット（スリット幅１２０μｍ、スペクトルの分解能１０ｃｍ−１に対応）上にラマン散乱光の像を結像し、刻線数６００ｇｒｏｏｖｅｓ／ｍｍの回折格子で分光して、電荷結合素子（Ｐｒｉｎｃｅｔｏｎ Ｉｎｓｔｒｕｍｅｎｔｓ，Ｉｎｃ．ＬＮ／ＣＣＤ−１０２４ＥＨＲＢ）にてラマンスペクトルを計測した。計測にあたり、図６Ａでは６００ｓの露光を１回、図６Ｈでは１００ｓの露光を６回行っている。なお、計測したデータの横軸は、電荷結合素子（１０２４個）のチャンネル番号であり、この状態でスパイクノイズを除去した後、横軸を波数に変換して図６Ａ〜６Ｉに示している。
図２乃至図５のアルゴリズム（本実施形態の方法）に相当するプログラムをソフトウェアＩｇｏｒ（ＷａｖｅＭｅｔｒｉｃｓ，Ｉｎｃ．）上で動作するよう作成し、ＲｅｃｏｇＰＡ，ＲａｔｉｏＰＡ，ＣＲＲを指定して、１０２４個のデータ列にスパイクノイズ除去処理を施した。図６Ａ（の横軸がチャンネル番号のデータ）に除去処理を施し、図６Ｂ（の横軸がチャンネル番号のデータ）を得るにあたり、計算に要する時間は約１０秒であった。なお、計算にはパーソナルコンピュータＩＢＭ Ｔｈｉｎｋ Ｐａｄ Ｘ２０を用いた。
データを計測するにあたり、前述の電荷結合素子の制御にソフトウェアＷｉｎｓｐｅｃ（Ｐｒｉｎｃｅｔｏｎ Ｉｎｓｔｒｕｍｅｎｔｓ，Ｉｎｃ．）を使用するが、必要に応じてスパイクノイズ除去機能を動作させることができる。図６Ａは、ソフトウェアＷｉｎｓｐｅｃにおけるスパイクノイズ除去機能を停止させ、計測したデータである。図６Ａでは、孤立したスパイクノイズあるいは２個以上連続したスパイクノイズが、多数観測されている。各種除去機能の比較のため、連続したスパイクノイズのうち２箇所に（＃），（＃＃）を施した。（＃）は全体の幅が８チャンネル（８４〜９１チャンネル）であり、幅３チャンネル、幅３チャンネル、幅４チャンネルの三つのスパイクノイズからなる。（＃＃）は全体の幅が８チャンネル（５５〜６２チャンネル）であり、幅４チャンネル、幅５チャンネルの二つのスパイクノイズからなる。なお、上記のスパイクノイズ幅の値は、図６Ａ〜６Ｈの横軸がチャンネル（波数に変換する前）の状態で読み取った結果である。
図６Ｂは、平均化反復法（本実施形態の方法）においてＲｅｃｏｇＰＡ＝８，ＲａｔｉｏＰＡ＝０．８，ＣＲＲ＝４とし、図６Ａのスパイクノイズ除去を試みた結果である。図６Ｂでは、波形が歪むことなく、全てのスパイクノイズが完全に除去されている。また図６Ｄは、本実施形態の方法（平均化反復法）においてＲｅｃｏｇＰＡ＝８，ＲａｔｉｏＰＡ＝０．８，ＣＲＲ＝１とし、図６Ａのスパイクノイズ除去を試みた結果である。図６Ｄでは波形の歪みはないものの、スパイクノイズが残存している（（＃），（＃＃）は、形状は変化するが、どちらも残存している）。ここで、閾値Ｙ１（ＳＩＧ×ＲｅｃｏｇＰＡ），Ｘ１−１（ＣＲＲ），Ｙ２（ＳＩＧ×ＲａｔｉｏＰＡ），Ｘ２（２×ＣＲＲ）は、ＲｅｃｏｇＰＡ，ＲａｔｉｏＰＡ，ＣＲＲで指定されるため、この３つのパラメータの値を示した。
ＣＲＲ＝４のとき、検出可能なスパイクノイズの幅は６チャンネルとなり（ステップＳ２２，Ｓ２６参照）、この値は（＃），（＃＃）を構成するスパイクノイズの幅より大きい。このため、（＃），（＃＃）を構成するスパイクノイズは、スパイクノイズとして検出され、除去処理の対象となる。その結果、（＃），（＃＃）も完全に除去されたと予想される。これに対しＣＲＲ＝１のとき、検出可能なスパイクノイズの幅は３チャンネルとなり、この値は（＃），（＃＃）を構成する一部のスパイクノイズの幅より小さい。このため、これらのスパイクノイズが検出されず、除去処理の対象にならなかったため、（＃），（＃＃）が残存したと考えられる。ＣＲＲ＝４のときとＣＲＲ＝１のときとの相違は、前述したＣＲＲの重要性を裏付けている。
ＣＲＲには、（スパイクノイズ幅）〜（スパイクノイズ幅）＋２を目安として指定するのが望ましい。（＃），（＃＃）を例にとれば、（＃），（＃＃）を構成するスパイクノイズの幅は最小３チャンネル、最大５チャンネルであるから、ＣＲＲは３〜７に指定することになる。ＲｅｃｏｇＰＡ，ＲａｔｉｏＰＡはＲｅｃｏｇＰＡ＝１０，ＲａｔｉｏＰＡ＝１を目安として指定するのが望ましい。本例では、ＣＲＲを３〜７の範囲で変化させ、それに応じてＲｅｃｏｇＰＡ，ＲａｔｉｏＰＡを調整した結果、連続したスパイクノイズを完全に除去し、かつデータ波形に歪みを発生しないような組み合わせとして、ＲｅｃｏｇＰＡ＝８，ＲａｔｉｏＰＡ＝０．８，ＣＲＲ＝４を見出した。
図６Ｅは、単純平均化法により、図６Ａのスパイクノイズ除去を試みた結果である。単純平均化法による計算を実現するにあたり、本実施形態の方法（平均化反復法）を示した図２において、ステップＳ１０，Ｓ２０の過程は変更せず、ステップＳ３０の過程のみを変更した。スパイクノイズ範囲ｊ≦ｉ≦ｊ＋ｑ＋１（ステップＳ２６参照）内のデータ値の全てを、スパイクノイズ範囲外のデータ値の平均値（ｙｊ−１＋ｙｊ＋ｑ＋２）／２で１回のみ置換することになる。単純平均化法の適用では、三つの閾値Ｙ１（ＳＩＧ×ＲｅｃｏｇＰＡ），Ｘ１−１（ＣＲＲ），Ｙ２（ＳＩＧ×ＲａｔｉｏＰＡ）を設定することに変更はなく、本実施形態の方法（平均化反復法）とスパイクノイズの検出法やスパイクノイズ範囲の決定法は同一である。そこでこれらの三つの閾値を図６Ｂの場合と同一にして、単純平均化法を適用した。図６Ｅから、波形の歪みはなくともスパイクノイズは残存していること、また（＃）は完全に除去されているが（＃＃）は残存していることがわかる。
図６Ｂ及び６Ｅを得るにあたり、ステップＳ１０，Ｓ２０は共通であるため、残存スパイクノイズの相違は、ステップＳ３０の変更の有無に起因する。図６Ｂ及び６Ｅの比較により、本実施形態の方法（平均化反復法）のスパイクノイズ除去過程において、平均化を１回で終了せず、複数のデータ値に対して平均化を反復することの意義が示された。
図６Ｆは、メディアンフィルター法を利用して、図６Ａのスパイクノイズ除去を試みた結果である。メディアンフィルター法に相当するプログラムをソフトウェアＩｇｏｒ（ＷａｖｅＭｅｔｒｉｃｓ，Ｉｎｃ．）上で動作するよう作成し、スパイクノイズ除去を試みた。ここでは幅７チャンネル内の値全てを、幅７チャンネル内の中央値で置換した。図６Ｆは、スパイクノイズが残存していること、及び（＃），（＃＃）の両者が残存していることを示している。また一部拡大図６Ｇの波形は、一部拡大図６Ｃに比較して大きく歪んでいることがわかる。
（＃），（＃＃）はいずれも幅８チャンネルのスパイクノイズであり、幅７チャンネル内の中央値を利用しても、これらを完全には除去できないことが、図６Ｆからわかった。図６Ｆにおける波形の歪みは、前述の予想を裏付けるものである。
図６Ｈは、ソフトウェアＷｉｎｓｐｅｃにおいて、時系列法を利用したスパイクノイズ除去機能を動作させながら、計測した結果である。時系列法の性質上、複数回の計測を行う必要がある。このため１００ｓの計測を６回行い、合計の計測時間を図１Ａの計測時間と同一の６００ｓとした。図６Ｈは、最終回の計測データであるが、このデータでは、１〜５回目の計測データの時間的相関性が考慮されている。図６Ｈは、スパイクノイズが完全に除去されたことを示している。一方で縦軸の値に注意し、一部拡大図６Ｉを一部拡大図６Ｃと比較すると、図６Ｉは図６Ｃに対して、信号雑音比が低下していることがわかる。
時系列法を適用したとき、スパイクノイズ除去効果は本実施形態の方法（平均化反復法）とほぼ同等であるが、本実施形態の方法（平均化反復法）に対し、信号雑音比が低下する。これは１回あたりの計測時間が短いためである。計測時間が短くなれば、混入するスパイクノイズの量は減少するが、シグナル強度も減少し、信号雑音比は低下する。
以上の詳説から、波形に歪みを発生させず、信号雑音比を確保して、連続したスパイクノイズを除去するにあたり、最も効果的な方法は、本実施形態の方法（平均化反復法）であることが示された。
本実施形態の方法（平均化反復法）は、第一にはスペクトルデータ等の２次元データ波形を処理するために考案されたが、これを、２次元の横軸に対して１価の縦軸を持つ３次元データ波形の処理に拡張することができる。すなわち、３次元データ波形に混入したスパイクノイズの除去に、本実施形態における１次微分値を１次偏微分値に変更してスパイクノイズ除去方法を適用することが可能である。その場合、ある１軸に沿った微分（偏微分）によってスパイクノイズ範囲を見出し、その近傍（２次元面内）で極値を探してスパイクノイズの中心位置を求める。その後、スパイクノイズ中心から距離１１／２の点、２１／２の点、４１／２の点、５１／２の点・・・とスパイクノイズの幅の２倍を超えない距離で、スパイクノイズ中心の周囲の縦軸の平均値を計算してゆき、ステップＳ３５，Ｓ３５ａと同様に収束の判定を行い、収束した際の平均値でデータを置換する。上記のような、平均値の計算と平均値による置換を、スパイクノイズ中心からスパイクノイズ幅以内の距離の点について反復する。さらに３次元以上の横軸に対して１価の縦軸を持つ４次元以上のデータ波形ついても、原理的にはこの方法を拡張することが可能である。
本実施形態は、スパイクノイズ検出・除去に必要な閾値の決定、スパイクノイズの検出とスパイクノイズ範囲の決定、スパイクノイズの除去という３つの過程（ステップ）からなる。
より詳細には、本実施形態のスパイクノイズ除去方法は、物理的測定により得られた各点の上のデータについて、データ波形（ｉ，ｙｉ）の（ａ）１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが正値でΔｙｇ，ｇ＋１が、又は（ｂ）Δｙｇ，ｇ＋１が負値で｜Δｙｇ，ｇ＋１｜が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、前記始点から所定の第１間隔Ｘ１の範囲内で、（ａ）１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、又は（ｂ）Δｙｈ，ｈ＋１が正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、（１）前記始点から前記終点までのスパイクノイズ範囲における（ａ）最大値（ｙｍ）又は（ｂ）最小値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、（２）（ａ）最大値（ｙｍ）又は（ｂ）最小値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、（３）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、（ａ）最大値（ｙｍ）又は（ｂ）最小値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、（４）データ値ｙｍ−ｃ，ｙｍ＋ｄ（ただしｃ，ｄは１≦ａ＜ｃ≦Ｘ２，１≦ｂ＜ｄ≦Ｘ２を満たす特定の値、Ｘ２は所定の第２間隔）の平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２を計算するステップと、（５）前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２と前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２との差を計算するステップと、（６）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２で置換するステップと、（７）文字ａの値を文字ｃの値で置換し、文字ｂの値を文字ｄの値で置換するステップと、（８）ａ＜ｃ≦Ｘ２，ｂ＜ｄ≦Ｘ２の範囲の各ｃ，各ｄに対して、前記（４）（５）（６）（７）を行い、前記スパイクノイズ点（ｍ，ｙｍ）に対する処理が終了したと判断するステップ（Ｓ３３）と、（９）前記スパイクノイズ点（ｍ，ｙｍ）から所定第３間隔Ｘ１−１の範囲内に位置する各データ点に対して、前記スパイクノイズ点（ｍ，ｙｍ）に対する、前記（１）（２）（３）（４）（５）（６）（７）（８）の処理と同様の処理を行い、当該スパイクノイズ範囲内のスパイクノイズが除去されたと判断するステップ（Ｓ３２）とを有する。
図１Ａ及び１Ｂは、スパイクノイズの除去の過程の実施形態の考え方を示す。
図１Ａは、スパイクノイズの除去の操作の説明のために作成した２次元データ波形であり、スパイクノイズ範囲を４≦ｉ≦７とし、ｙ５＝５０，ｙ６＝４６を、バックグラウンドに混入した２本の連続したスパイクノイズと仮定する。図中折線（ａ１）は生データ、折線（ａ２）はスパイクノイズ除去後のデータを示す。
また、前記閾値Ｘ１−１を３、計算終了判定のための閾値（前記閾値Ｙ２）として、ＮＬ（ノイズレベル）＝４と仮定し、図１Ｂには、スパイクノイズを含む生データが、スパイクノイズの除去過程で変化してゆく様子を示した。
図１Ｂ第２列において、生データのスパイクノイズ範囲４≦ｉ≦７内の最大値はｙ５＝５０であり、これが最初の除去対象となりうる。ｙ５＝５０と平均値ｙ５１＝（ｙ４＋ｙ６）／２＝２５との差を計算すると｜ｙ５−ｙ５１｜＝２５となり、この値は（ノイズレベル）＝４よりも大きい。そこで、ｙ５＝５０をｙ５１＝２５で置換する。次に、ｙ５１＝２５と平均値ｙ５２＝（ｙ３＋ｙ７）／２＝２との差を計算すると｜ｙ５２−ｙ５１｜＝２３となり、この値は（ノイズレベル）＝４よりも大きい。そこで、ｙ５１＝２５をｙ５２＝２で置換する。さらに、ｙ５２＝２と平均値ｙ５３＝（ｙ２＋ｙ８）／２＝３との差を計算すると｜ｙ５３−ｙ５２｜＝１となり、この値は（ノイズレベル）＝４よりも小さい。そこで、スパイクノイズｙ５＝５０の除去が終了したと判断し、ｙ５２＝２で置換した状態で計算を終了する。この段階でのデータ波形は、図１Ｂ第３列に示した通りである。上記と同様の処理を、ｉ＝５，５−１，５＋１，５−２，５＋２，５−３，５＋３＝５，４，６，３，７，２，８なるｙｉについて行う（図５のステップＳ２２，Ｓ３２において、Ｘ１＝４とした場合に相当する）。ｙ５に関する処理は終了したので、ｙ５に隣接するｙ４に関する処理を説明する。
ｙ４＝５とｙ４１＝（ｙ３＋ｙ５）／２＝２．５との差を計算すると、｜ｙ４−ｙ４１｜＝２．５となり、この値は（ノイズレベル）＝４よりも小さい。そこで、ｙ４＝５をスパイクノイズとはみなさずに、ｙ４＝５のまま計算を終了する。これにより、データ波形は図１Ｂ第３列の状態を保持し、図１Ｂ第４列のようになる。次に、ｙ４と同様にｙ５に隣接するｙ６に関する処理を説明する。
ｙ６＝４５と平均値ｙ６１＝（ｙ５＋ｙ７）／２＝１．５との差を計算すると｜ｙ６−ｙ６１｜＝４３．５となり、この値は（ノイズレベル）＝４よりも大きい。そこで、ｙ６＝４５をｙ６１＝１．５で置換する。次に、ｙ６１＝１．５と平均値ｙ６２＝（ｙ４＋ｙ８）／２＝５との差を計算すると｜ｙ６２−ｙ６１｜＝３．５となり、この値は（ノイズレベル）＝４よりも小さい。そこで、スパイクノイズｙ６＝４５の除去が終了したと判断し、ｙ６１＝１．５で置換した状態で計算を終了する。これにより、データ波形は、図１Ｂ第５列のように変化する。
上記ｙ５，ｙ４，ｙ６に関する処理と同様の処理を、ｙ３，ｙ７，ｙ２，ｙ８にも施すと、ｙ４の場合と同様に個々の値は変化せず、最終的なデータ波形は図１Ｂ第６列のようになる。以上の過程により、生データにおいて連続する２本のスパイクノイズｙ５＝５０，ｙ６＝４５が完全に除去され、平均値による置換を、範囲を拡大しながら反復することの利点が示された。
なお、単純平均化法（特許第２６２９５１１号の方法）では、生データは図１Ｂ第７列のように変化する。特許第２６２９５１１号の方法では、スパイクノイズ範囲は５≦ｉ≦６となり、スパイクノイズ範囲を挟むｙ４＝５，ｙ７＝１の平均値３でスパイクノイズ範囲内の値が全て置換されるためである。このように単純平均化法では、一部のデータを全て同一の値で置換するため、波形に歪みを発生することが予想される。また、メディアンフィルター法で３チャンネル内の中央値による置換を仮定すると、生データのｙ０，ｙ１，ｙ２，ｙ３，ｙ４，ｙ５，ｙ６，ｙ７，ｙ８は順に０，１，３，３，５，４５，４５，５，４，４，２に変化し、スパイクノイズが残存することになる。単純平均化法で発生が予想される波形の歪みは、図１のようにバックグラウンド上ではなく、シグナル波形上にスパイクノイズが混入した場合に顕著である。具体例として、特許第２６２９５１１号図８の測定時間３．２秒付近の波形に観測されるショルダーピークが挙げられる。本実施形態では、スパイクノイズ近傍のシグナルデータ値の平均値で置換を行う。このため置換後の波形には、スパイクノイズ近傍のシグナルデータの特徴が反映され、単純平均化法、メディアンフィルター法の場合よりも、置換後の波形の歪みが小さくなる。
図２〜図５は、本発明の実施形態を示すフローチャートである。図２に示すように、ステップＳ１０でスパイクノイズの検出・除去で必要となる閾値の決定を行い、ステップＳ２０でスパイクノイズの検出とスパイクノイズ範囲の決定を行い、ステップＳ３０はスパイクノイズの除去を行う。
後述するように、オペレータは、パラメータとして、ステップＳ１０で、バックグラウンド上の２点の位置、および認識パラメータ（ＲｅｃｏｇＰＡ）、およびノイズ比パラメータ（ＲａｔｉｏＰＡ）を決定し、Ｓ２０でスパイクノイズレンジ（ＣＲＲ）を決定する。
バックグラウンド上の上記の２点の回帰直線とデータの差の標準偏差がＳＩＧとして算出され、前記閾値Ｙ１に相当するのはＮＴ＝ＳＩＧ×ＲｅｃｏｇＰＡとなり、前記閾値Ｘ１−１に相当するのはＣＲＲとなり、前記閾値Ｙ２に相当するのはＮＬ＝ＳＩＧ×ＲａｔｉｏＰＡとなる。なお、本発明の実施形態はこのような閾値の決定法に限定されるものではなく、別の方法により前記閾値Ｙ１およびＹ２を指定することも可能である。
ステップＳ１０、Ｓ２０でパラメータが、一度指定されると、該パラメータはステップＳ２０、Ｓ３０で利用され、計算は自動的に進行する。
なお、図２〜図５において、オペレータによる指定が含まれる過程は平行四辺形で表示され、それ以外の過程は長方形あるいはひし形で表示されている。なお、ひし形は条件による分岐を意味している。また、連続するデータ間での１次微分値を差分と定義し、差分を用いて計算過程を説明してある。
以下、図３〜図５のフローチャートを図１のスパイクノイズ除去過程を参照して説明する。図３のステップＳ１１〜Ｓ１６は図２のステップＳ１０を詳細に示し、図４のステップＳ２１〜Ｓ２６、Ｓ２５ａは図２のステップＳ２０を詳細に示し、図５のステップＳ３１〜Ｓ３５、Ｓ３５ａは図２のステップＳ３０を詳細に示す。
ここで、前記閾値をＹ１＝１２（ＳＩＧ×ＲｅｃｏｇＰＡ＝１２），Ｘ１＝４（ＣＲＲ＝３），Ｙ２＝４（ＳＩＧ×ＲａｔｉｏＰＡ＝４、図１Ａにおけるノイズレベル），Ｘ２＝６（２×ＣＲＲ＝６）と仮定する。後述するように、これらの値は、図３のステップＳ１５、Ｓ１６、図４のステップＳ２２で指定される。
まず、図３に示すスパイクノイズの検出・除去で必要となる閾値の決定について説明する。
ステップＳ１１では、バックグラウンド上の２点を指定し、ステップＳ１２では、２点間の回帰直線とデータ値の差の標準偏差を計算してＳＩＧとする。ステップＳ１３で認識パラメータＲｅｃｏｇＰＡを指定し、ステップＳ１４でノイズ比パラメータＲａｔｉｏＰＡを指定する。ステップＳ１５では、ＮＴ＝ＳＩＧ×ＲｅｃｏｇＰＡを閾値（Ｙ１）として算出し、ステップＳ１６では、ＮＴ＝ＳＩＧ×ＲａｔｉｏＰＡを閾値（Ｙ２）として算出する。
次に、図４に示すスパイクノイズの検出とスパイクノイズ範囲の決定について説明する。
ステップＳ２１では、検出するスパイクノイズピーク値の正負を指定する（以下、スパイクノイズピーク値が正の場合を仮定する）。ステップＳ２２では、検出・除去の対象として許容するスパイクノイズ幅のチャンネル数の最大値より、２減じた値を、スパイクノイズレンジ（ＣＲＲ）として指定する。ここで、ＣＲＲ＝Ｘ１−１とする。ステップＳ２３では、データ（ｉ，ｙｉ）の差分Δｙｉ，ｉ＋１＝ｙｉ＋１−ｙｉを計算する。
パラメータの指定および計算がステップＳ２３まで進行した後、ステップＳ２４において、Δｙｊ，ｊ＋１＝ｙｊ＋１−ｙｊ＞＋ＮＴのとき、ｊをスパイクノイズの始まるチャンネルとする。図１Ａ及び１Ｂの例の場合、Δｙ４，５＝ｙ５−ｙ４＝４５＞１２であるのでｊ＝４である。
ステップＳ２５においてΔｙｊ，ｊ＋ｑ＋１＝ｙｊ＋ｑ＋１−ｙｊ＜−ＮＴのとき、ｊ＋ｑ＋１をスパイクノイズが終了するチャンネルとする。ただし、１≦ｑ≦ＣＲＲとする。図１の例の場合、Δｙ６，７＝ｙ７−ｙ６＝−４４＜−１２であり、ｊ＋ｑ＋１＝７であるから、ｑ＝２となる。
ステップＳ２６において、スパイクノイズ範囲はｊ≦ｉ≦ｊ＋ｑ＋１であるから、スパイクノイズ幅はチャンネル数ｑ＋２（間隔ｑ＋１）である。ただし、１≦ｑ≦ＣＲＲである。図１の例の場合、スパイクノイズ範囲は４≦ｉ≦７、ｑ＝２となる。ここでＣＲＲ＝３であるから、１≦ｑ≦ＣＲＲが成立する。ここで、スパイクノイズ範囲は４〜７チャンネルであり、スパイクノイズ幅は４チャンネルである。
ステップＳ２５ａでは、１≦ｑ≦ＣＲＲの範囲で差分がＳ２５の条件を満たさないときは、スパイクノイズは存在しないとみなし、計算を終了する。このステップＳ２５ａの過程は、スパイクノイズ幅が過剰に大きくなることを防止するために設けた。
ステップＳ２５ａの過程を省略すると、スパイクノイズ幅が、正常な信号ピークのバンド幅に相当する値を持つ場合、正常な信号ピークを広範囲におよぶスパイクノイズと判断し、このピークを図５のステップＳ３１以下の過程で除去してしまう可能性がある。
図１Ａ及び１Ｂの例の場合、ＣＲＲ＝３であるからスパイクノイズ幅のチャンネル数ｑ＋２は５以下となる必要がある。実際にはｑ＋２＝４であり、上記制限内であることがわかる。このようにＣＲＲの設定は、スパイクノイズ幅を制限し、連続したスパイクノイズと、正常な信号ピークとを区別する機能を持つ。
次に、図５に示すスパイクノイズの除去について説明する。
ステップＳ３１において、スパイクノイズ範囲ｊ≦ｉ≦ｊ＋ｑ＋１内でｙｉが最大の点（ｐ，ｙｐ）を特定する。図１の例の場合、スパイクノイズ範囲４≦ｉ≦７内でｙｉが最大の点は（５，ｙ５）であるので、ｐ＝５である。
ステップＳ３２では、ステップＳ３３〜Ｓ３５、Ｓ３５ａをｉ＝ｐ，ｐ−１，ｐ＋１，‥，ｐ−ＣＲＲ，ｐ＋ＣＲＲについて実行する。
ステップＳ３３では、ステップＳ３４、Ｓ３５、Ｓ３５ａを０≦ｒ≦２×ＣＲＲについて実行する。ここで、２×ＣＲＲ＝Ｘ２である。図１の例の場合、前記（５，ｙ５）を処理するにあたり、ｒは２×ＣＲＲ＝６以下となる必要がある。本例では、前述の通り、ｒ＝０，１，２で（５，ｙ５）の処理が完了しており、これは上記制限内であることがわかる。
ステップＳ３４では、ｙｉｒ＝（ｙｉ−ｒ＋ｙｉ＋ｒ）／２を計算する。
ステップＳ３４の後処理は分岐し、｜ｙｉｒ＋１−ｙｉｒ｜＞ＮＬの場合に処理はステップＳ３５に進み、｜ｙｉｒ＋１−ｙｉｒ｜＜ＮＬの場合に処理はステップＳ３５ａに進む。
このようにして（５，ｙ５）の処理がステップＳ３５ａで完了した後、再度ステップＳ３２に戻ることになる。図１Ａ及び１Ｂの例の場合、ステップＳ３２〜Ｓ３５ａで、上記と同様の処理を、ｉ＝５，５−１，５＋１，５−２，５＋２，５−３，５＋３＝５，４，６，３，７，２，８なるｙｉに施すことになる。そこで、ｙ５に関する説明の次に、ｙ５に隣接するｙ４に関する説明を行う。
ステップＳ３３〜Ｓ３５ａで、（４，ｙ４）を処理するにあたり、ｒは２×ＣＲＲ＝６以下となる必要がある。本例では、前述の通り、ｒ＝０で（４，ｙ４）の処理が完了しており、これは上記制限内であることがわかる。このようにして（４，ｙ４）の処理がステップＳ３５ａで完了し、再度ステップＳ３２に戻ることになる。次にｙ４と同様にｙ５に隣接するｙ６の処理に関する説明を行う。
ステップＳ３３〜Ｓ３５ａで、（６，ｙ６）を処理するにあたり、ｒは２×ＣＲＲ＝６以下となる必要がある。本例では、前述の通り、ｒ＝０，１で（６，ｙ６）の処理が完了しており、これは上記制限内であることがわかる。このようにして（６，ｙ６）の処理がステップＳ３５ａで完了した後、再度ステップＳ３２に戻ることになる。
ｙ５，ｙ４，ｙ６に関する処理と同様の処理をｙ３，ｙ７，ｙ２，ｙ８にも施すと、前述の通り、ｙ４の場合と同様に個々の値は変化せずｒ＝０で処理は完了し、ｒが２×ＣＲＲ＝６以下という制限内に収まっていることがわかる。このようにして、図１に示したように、連続した２本のスパイクノイズが除去されることになる。
本実施形態の本質は、平均値による置換の反復にあるが、これには二つの意味がある。第一は、図５のステップＳ３２の過程で、スパイクノイズ範囲内の最大値のほか、最大値近傍の複数の点が除去処理対象となること、第二は、図５のステップＳ３３の過程で、除去処理対象となった個々の点の値が、当該点周辺のベースラインに近似できるまで置換を繰り返すことである。
前記閾値Ｙ１，Ｘ１−１，Ｙ２，Ｘ２は、これまでの記述通り、ＳＩＧ×ＲｅｃｏｇＰＡ，ＣＲＲ，ＳＩＧ×ＲａｔｉｏＰＡ，２×ＣＲＲとするほか、別の方法でも指定できる。一例として、Ｘ２を１．５×ＣＲＲ，３×ＣＲＲ，４×ＣＲＲなどに設定してもよい。前記閾値Ｙ１，Ｘ１−１，Ｙ２，Ｘ２の指定法、あるいは上記のＸ２の指定法では、Ｘ２≧Ｘ１−１となる。特定のデータ値を平均値の置換で小さくする（ステップＳ３３）にあたり、Ｘ２≧Ｘ１−１は、重要な意味を持つ。Ｘ２≧Ｘ１−１の持つ意味は、上記の特定のデータ値から許容されるスパイクノイズの幅のチャンネル数（より２を減じた値）以上の距離内のデータを用いて、平均値を算出することにある。許容されるスパイクノイズの幅のチャンネル数（より２を減じた値）を下回る距離内のデータのみを用い、平均値を算出しても、スパイクノイズ範囲内のデータ値の平均値を利用することになり、スパイクノイズの除去が不十分になる可能性が高い。従って、前述の方法とは異なる方法で閾値Ｙ１，Ｘ１−１，Ｙ２，Ｘ２を指定する場合も、Ｘ２≧Ｘ１−１となることが望ましい。
閾値Ｘ１−１（ＣＲＲ）は、ステップＳ２２の通り、スパイクノイズ幅のチャンネル数の最大値より２減じた値を意味する。スパイクノイズ幅のチャンネル数の最大値がＸ１＋１となるから、スパイクノイズ範囲の始点から終点に至る間隔の最大値は、Ｘ１となる。Ｘ１は、請求項１記載の第１間隔に相当する。また、閾値Ｘ２（２×ＣＲＲ）は、ステップＳ３３の通り、除去処理対象となる特定のスパイクノイズのチャンネル番号と、平均値の算出に使用するデータのチャンネル番号との差の最大値を意味する。前記のデータのチャンネルと、前記のスパイクノイズのチャンネルとの差の最大値は、間隔にしてＸ２となる。Ｘ２は、請求項１、請求項２記載の第２間隔に相当する。さらにステップＳ３２の通り、除去処理対象となるデータ値のチャンネル番号と、スパイクノイズ範囲内の最大値のチャンネル番号との差の最大値は、Ｘ１−１（ＣＲＲ）である。これは、間隔にしてＸ１−１となり、請求項３記載の第３間隔に相当する。
前述の方法で閾値Ｙ１（ＳＩＧ×ＲｅｃｏｇＰＡ），Ｘ１−１（ＣＲＲ），Ｙ２（ＳＩＧ×ＲａｔｉｏＰＡ），Ｘ２（２×ＣＲＲ）を指定するにあたり、除去処理後も連続したスパイクノイズが残存するときはＸ１−１（ＣＲＲ）の値を大きくし、データ波形が歪むときはＸ１−１（ＣＲＲ）の値を小さくすればよい。閾値Ｘ１−１（ＣＲＲ）の値を大きくすること、ステップＳ３２が示すように、より多数のデータ値が除去処理の対象となる。また同時に閾値Ｘ２（２×ＣＲＲ）の値も大きくなり、ステップＳ３３が示すように、平均値の算出でより広範囲のデータ値を利用できるので、除去処理対象となった特定のデータ値を、スパイクノイズ範囲のデータ値の影響を受けずに、当該点周辺のベースラインに近似した値に置換できる。このように、閾値Ｘ１−１（ＣＲＲ）の値を大きくすることで、スパイクノイズのうち、高い連続性を示すもののを効果的に除去できる。上記の影響を考慮すれば、閾値Ｘ１−１（ＣＲＲ）は四つの閾値の中で最も本質的である。一方でステップＳ３２に代わり、ｉ＝ｐ，ｐ−１，ｐ＋１，…，ｐ−Ｘ３，ｐ＋Ｘ３とし、スパイクノイズの連続性が高いとき、Ｘ１−１，Ｘ２とは独立にＸ３を大きな値に設定して計算を実施してもよい。
また、滑らかなバックグラウンド上のスパイクノイズを除去するにあたり、除去処理後もスパイクノイズが残存するときは、ＲｅｃｏｇＰＡの値を大きくするのがよい。その理由は、滑らかなバックグラウンド上では標準偏差ＳＩＧの値が小さいために閾値Ｙ１（ＳＩＧ×ＲｅｃｏｇＰＡ）の値が小さくなり、ステップＳ２４、Ｓ２５でのスパイクノイズ範囲の決定が困難となるためである。
前述の方法で閾値Ｙ１，Ｘ１−１，Ｙ２，Ｘ２を指定するにあたり、連続したスパイクノイズが、滑らかなバックグラウンド上に存在する場合も含めて完全に除去され、かつデータ波形に歪みが発生しないよう各閾値を指定する必要があるが、前述の通り、前記の指定は過度の実験を必要としない。
以下、図６Ａ〜６Ｉに示した実施例により本実施形態を説明する。図６Ａは、ポリメタクリル酸メチルからなる板（以下、アクリル板と称する）のラマンスペクトルを、電荷結合素子を用い、スパイクノイズ除去機能なし、計測時間６００秒で計測した結果であり、宇宙線スパイクノイズが混入したスペクトルとなっている。図６Ｂ，６Ｄ，６Ｅ，６Ｆは、図６Ａに対して各種のスパイクノイズ除去方法を適用した結果である。図６Ｈは測定と同時に、時系列法によるスパイクノイズ除去機能を動作させ、スパイクノイズ除去を試みた結果である。図６Ｃ，６Ｇ，６Ｉは、それぞれ図６Ｂ，６Ｆ，６Ｈの一部拡大図である。ラマンスペクトルは、試料にレーザー光を照射し、試料が発する散乱光のうち照射光よりも低エネルギーのものについて、エネルギーの損失量毎に、散乱光の強度をプロットした２次元データである。エネルギーの損失量は、化学結合の種類や状態に依存するため、ここで得られたラマンスペクトルは、アクリル板固有のスペクトルとなる。
図６Ａ及び６Ｈのデータは、ダイオード励起可視レーザー（Ｓｐｅｃｔｒａ Ｐｈｙｓｉｃｓ Ｍｉｌｌｅｎｉａ）から発振する波長５３２ｎｍのレーザー光（強度０．６ｍＷ，ビーム径約０．２ｍｍ）を、厚さ２ｍｍのアクリル板に、アクリル板法線から約４５°の角度で照射して測定した。アクリル板から発生するラマン散乱光とレイリー散乱光を、アクリル板法線方向に設置したＦ／１の集光レンズで集光して平行光とした後、ノッチフィルター（Ｋａｉｓｅｒ Ｈｏｌｏｇｒａｐｈｉｃ Ｎｏｔｃｈ Ｆｉｌｔｅｒ ＨＮＰＦ−５３２．０−１．５）を通過させてレイリー散乱光をカットした。その後、Ｆ／７．５の集光レンズでシングルポリクロメータ（Ｊｏｂｉｎ Ｙｖｏｎ Ｒａｍａｎｏｒ Ｔ−６４０００，Ｆ／７．５）のスリット（スリット幅１２０μｍ、スペクトルの分解能１０ｃｍ−１に対応）上にラマン散乱光の像を結像し、刻線数６００ｇｒｏｏｖｅｓ／ｍｍの回折格子で分光して、電荷結合素子（Ｐｒｉｎｃｅｔｏｎ Ｉｎｓｔｒｕｍｅｎｔｓ，Ｉｎｃ．ＬＮ／ＣＣＤ−１０２４ＥＨＲＢ）にてラマンスペクトルを計測した。計測にあたり、図６Ａでは６００ｓの露光を１回、図６Ｈでは１００ｓの露光を６回行っている。なお、計測したデータの横軸は、電荷結合素子（１０２４個）のチャンネル番号であり、この状態でスパイクノイズを除去した後、横軸を波数に変換して図６Ａ〜６Ｉに示している。
図２乃至図５のアルゴリズム（本実施形態の方法）に相当するプログラムをソフトウェアＩｇｏｒ（ＷａｖｅＭｅｔｒｉｃｓ，Ｉｎｃ．）上で動作するよう作成し、ＲｅｃｏｇＰＡ，ＲａｔｉｏＰＡ，ＣＲＲを指定して、１０２４個のデータ列にスパイクノイズ除去処理を施した。図６Ａ（の横軸がチャンネル番号のデータ）に除去処理を施し、図６Ｂ（の横軸がチャンネル番号のデータ）を得るにあたり、計算に要する時間は約１０秒であった。なお、計算にはパーソナルコンピュータＩＢＭ Ｔｈｉｎｋ Ｐａｄ Ｘ２０を用いた。
データを計測するにあたり、前述の電荷結合素子の制御にソフトウェアＷｉｎｓｐｅｃ（Ｐｒｉｎｃｅｔｏｎ Ｉｎｓｔｒｕｍｅｎｔｓ，Ｉｎｃ．）を使用するが、必要に応じてスパイクノイズ除去機能を動作させることができる。図６Ａは、ソフトウェアＷｉｎｓｐｅｃにおけるスパイクノイズ除去機能を停止させ、計測したデータである。図６Ａでは、孤立したスパイクノイズあるいは２個以上連続したスパイクノイズが、多数観測されている。各種除去機能の比較のため、連続したスパイクノイズのうち２箇所に（＃），（＃＃）を施した。（＃）は全体の幅が８チャンネル（８４〜９１チャンネル）であり、幅３チャンネル、幅３チャンネル、幅４チャンネルの三つのスパイクノイズからなる。（＃＃）は全体の幅が８チャンネル（５５〜６２チャンネル）であり、幅４チャンネル、幅５チャンネルの二つのスパイクノイズからなる。なお、上記のスパイクノイズ幅の値は、図６Ａ〜６Ｈの横軸がチャンネル（波数に変換する前）の状態で読み取った結果である。
図６Ｂは、平均化反復法（本実施形態の方法）においてＲｅｃｏｇＰＡ＝８，ＲａｔｉｏＰＡ＝０．８，ＣＲＲ＝４とし、図６Ａのスパイクノイズ除去を試みた結果である。図６Ｂでは、波形が歪むことなく、全てのスパイクノイズが完全に除去されている。また図６Ｄは、本実施形態の方法（平均化反復法）においてＲｅｃｏｇＰＡ＝８，ＲａｔｉｏＰＡ＝０．８，ＣＲＲ＝１とし、図６Ａのスパイクノイズ除去を試みた結果である。図６Ｄでは波形の歪みはないものの、スパイクノイズが残存している（（＃），（＃＃）は、形状は変化するが、どちらも残存している）。ここで、閾値Ｙ１（ＳＩＧ×ＲｅｃｏｇＰＡ），Ｘ１−１（ＣＲＲ），Ｙ２（ＳＩＧ×ＲａｔｉｏＰＡ），Ｘ２（２×ＣＲＲ）は、ＲｅｃｏｇＰＡ，ＲａｔｉｏＰＡ，ＣＲＲで指定されるため、この３つのパラメータの値を示した。
ＣＲＲ＝４のとき、検出可能なスパイクノイズの幅は６チャンネルとなり（ステップＳ２２，Ｓ２６参照）、この値は（＃），（＃＃）を構成するスパイクノイズの幅より大きい。このため、（＃），（＃＃）を構成するスパイクノイズは、スパイクノイズとして検出され、除去処理の対象となる。その結果、（＃），（＃＃）も完全に除去されたと予想される。これに対しＣＲＲ＝１のとき、検出可能なスパイクノイズの幅は３チャンネルとなり、この値は（＃），（＃＃）を構成する一部のスパイクノイズの幅より小さい。このため、これらのスパイクノイズが検出されず、除去処理の対象にならなかったため、（＃），（＃＃）が残存したと考えられる。ＣＲＲ＝４のときとＣＲＲ＝１のときとの相違は、前述したＣＲＲの重要性を裏付けている。
ＣＲＲには、（スパイクノイズ幅）〜（スパイクノイズ幅）＋２を目安として指定するのが望ましい。（＃），（＃＃）を例にとれば、（＃），（＃＃）を構成するスパイクノイズの幅は最小３チャンネル、最大５チャンネルであるから、ＣＲＲは３〜７に指定することになる。ＲｅｃｏｇＰＡ，ＲａｔｉｏＰＡはＲｅｃｏｇＰＡ＝１０，ＲａｔｉｏＰＡ＝１を目安として指定するのが望ましい。本例では、ＣＲＲを３〜７の範囲で変化させ、それに応じてＲｅｃｏｇＰＡ，ＲａｔｉｏＰＡを調整した結果、連続したスパイクノイズを完全に除去し、かつデータ波形に歪みを発生しないような組み合わせとして、ＲｅｃｏｇＰＡ＝８，ＲａｔｉｏＰＡ＝０．８，ＣＲＲ＝４を見出した。
図６Ｅは、単純平均化法により、図６Ａのスパイクノイズ除去を試みた結果である。単純平均化法による計算を実現するにあたり、本実施形態の方法（平均化反復法）を示した図２において、ステップＳ１０，Ｓ２０の過程は変更せず、ステップＳ３０の過程のみを変更した。スパイクノイズ範囲ｊ≦ｉ≦ｊ＋ｑ＋１（ステップＳ２６参照）内のデータ値の全てを、スパイクノイズ範囲外のデータ値の平均値（ｙｊ−１＋ｙｊ＋ｑ＋２）／２で１回のみ置換することになる。単純平均化法の適用では、三つの閾値Ｙ１（ＳＩＧ×ＲｅｃｏｇＰＡ），Ｘ１−１（ＣＲＲ），Ｙ２（ＳＩＧ×ＲａｔｉｏＰＡ）を設定することに変更はなく、本実施形態の方法（平均化反復法）とスパイクノイズの検出法やスパイクノイズ範囲の決定法は同一である。そこでこれらの三つの閾値を図６Ｂの場合と同一にして、単純平均化法を適用した。図６Ｅから、波形の歪みはなくともスパイクノイズは残存していること、また（＃）は完全に除去されているが（＃＃）は残存していることがわかる。
図６Ｂ及び６Ｅを得るにあたり、ステップＳ１０，Ｓ２０は共通であるため、残存スパイクノイズの相違は、ステップＳ３０の変更の有無に起因する。図６Ｂ及び６Ｅの比較により、本実施形態の方法（平均化反復法）のスパイクノイズ除去過程において、平均化を１回で終了せず、複数のデータ値に対して平均化を反復することの意義が示された。
図６Ｆは、メディアンフィルター法を利用して、図６Ａのスパイクノイズ除去を試みた結果である。メディアンフィルター法に相当するプログラムをソフトウェアＩｇｏｒ（ＷａｖｅＭｅｔｒｉｃｓ，Ｉｎｃ．）上で動作するよう作成し、スパイクノイズ除去を試みた。ここでは幅７チャンネル内の値全てを、幅７チャンネル内の中央値で置換した。図６Ｆは、スパイクノイズが残存していること、及び（＃），（＃＃）の両者が残存していることを示している。また一部拡大図６Ｇの波形は、一部拡大図６Ｃに比較して大きく歪んでいることがわかる。
（＃），（＃＃）はいずれも幅８チャンネルのスパイクノイズであり、幅７チャンネル内の中央値を利用しても、これらを完全には除去できないことが、図６Ｆからわかった。図６Ｆにおける波形の歪みは、前述の予想を裏付けるものである。
図６Ｈは、ソフトウェアＷｉｎｓｐｅｃにおいて、時系列法を利用したスパイクノイズ除去機能を動作させながら、計測した結果である。時系列法の性質上、複数回の計測を行う必要がある。このため１００ｓの計測を６回行い、合計の計測時間を図１Ａの計測時間と同一の６００ｓとした。図６Ｈは、最終回の計測データであるが、このデータでは、１〜５回目の計測データの時間的相関性が考慮されている。図６Ｈは、スパイクノイズが完全に除去されたことを示している。一方で縦軸の値に注意し、一部拡大図６Ｉを一部拡大図６Ｃと比較すると、図６Ｉは図６Ｃに対して、信号雑音比が低下していることがわかる。
時系列法を適用したとき、スパイクノイズ除去効果は本実施形態の方法（平均化反復法）とほぼ同等であるが、本実施形態の方法（平均化反復法）に対し、信号雑音比が低下する。これは１回あたりの計測時間が短いためである。計測時間が短くなれば、混入するスパイクノイズの量は減少するが、シグナル強度も減少し、信号雑音比は低下する。
以上の詳説から、波形に歪みを発生させず、信号雑音比を確保して、連続したスパイクノイズを除去するにあたり、最も効果的な方法は、本実施形態の方法（平均化反復法）であることが示された。
本実施形態の方法（平均化反復法）は、第一にはスペクトルデータ等の２次元データ波形を処理するために考案されたが、これを、２次元の横軸に対して１価の縦軸を持つ３次元データ波形の処理に拡張することができる。すなわち、３次元データ波形に混入したスパイクノイズの除去に、本実施形態における１次微分値を１次偏微分値に変更してスパイクノイズ除去方法を適用することが可能である。その場合、ある１軸に沿った微分（偏微分）によってスパイクノイズ範囲を見出し、その近傍（２次元面内）で極値を探してスパイクノイズの中心位置を求める。その後、スパイクノイズ中心から距離１１／２の点、２１／２の点、４１／２の点、５１／２の点・・・とスパイクノイズの幅の２倍を超えない距離で、スパイクノイズ中心の周囲の縦軸の平均値を計算してゆき、ステップＳ３５，Ｓ３５ａと同様に収束の判定を行い、収束した際の平均値でデータを置換する。上記のような、平均値の計算と平均値による置換を、スパイクノイズ中心からスパイクノイズ幅以内の距離の点について反復する。さらに３次元以上の横軸に対して１価の縦軸を持つ４次元以上のデータ波形ついても、原理的にはこの方法を拡張することが可能である。
Claims (9)
- 物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが正値でΔｙｇ，ｇ＋１が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、
前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、
（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最大値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、
（２）前記最大値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、
（３）前記差が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差が、当該所定閾値以上であるとき、最大値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、
を有するスパイクノイズ除去方法。 - 物理的測定により得られた各点の上のデータｙｉ（ｉは０以上の整数）について、データ波形（ｉ，ｙｉ）の１次差分値又は１次微分値Δｙｇ，ｇ＋１＝ｙｇ＋１−ｙｇが負値で｜Δｙｇ，ｇ＋１｜が、所定の正の第１閾値（Ｙ１）以上の値となるとき、データ点（ｇ，ｙｇ）をスパイクノイズ範囲の始点と判断するステップ（Ｓ２４）と、
前記始点から所定の第１間隔Ｘ１の範囲内で、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上の値となるとき、データ点（ｈ＋１，ｙｈ＋１）をスパイクノイズ範囲の終点と判断するステップ（Ｓ２５）と、
（１）前記始点から前記終点までのスパイクノイズ範囲における前記データ値の最小値（ｙｍ）を挟むデータ値ｙｍ−ａ，ｙｍ＋ｂ（ただしａ，ｂは１≦ａ，ｂ≦Ｘ２を満たす特定の値、Ｘ２は所定第２間隔）の平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を計算するステップと、
（２）前記最小値（ｙｍ）と前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２との差を計算するステップと、
（３）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、最小値（ｙｍ）を平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２で置換するステップと、
を有するスパイクノイズ除去方法。 - ステップ（３）で前記差の絶対値が、当該所定閾値以上であるとき、前記ステップ（１）（２）（３）の後に、
（４）データ値ｙｍ−ｃ，ｙｍ＋ｄ（ただしｃ，ｄは１≦ａ＜ｃ≦Ｘ２，１≦ｂ＜ｄ≦Ｘ２を満たす特定の値、Ｘ２は所定の第２間隔）の平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２を計算するステップと、
（５）前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２と前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２との差を計算するステップと、
（６）前記差の絶対値が、所定の第２閾値（Ｙ２）以下であるとき、処理を終了し、前記差の絶対値が、当該所定閾値以上であるとき、前記平均値（ｙｍ−ａ＋ｙｍ＋ｂ）／２を前記平均値（ｙｍ−ｃ＋ｙｍ＋ｄ）／２で置換するステップと、
（７）値ａを値ｃで置換し、値ｂを値ｄで置換するステップと、
（８）ａ＜ｃ≦Ｘ２，ｂ＜ｄ≦Ｘ２の範囲の各ｃ，各ｄに対して、前記（４）（５）（６）（７）を行うステップ（Ｓ３３）と、
をさらに有する請求の範囲第１項又は第２項に記載のスパイクノイズ除去方法。 - 前記ステップ（１）（２）（３）（４）（５）（６）（７）（８）の後に、
（９）前記スパイクノイズ点（ｍ，ｙｍ）から所定第３間隔Ｘ１−１の範囲内に位置する各データ点に対して、前記スパイクノイズ点（ｍ，ｙｍ）に対する、前記（１）（２）（３）（４）（５）（６）（７）（８）の処理と同様の処理を行い、当該スパイクノイズ範囲内のスパイクノイズが除去されたと判断するステップ（Ｓ３２）と、
をさらに有する請求の範囲第３項に記載のスパイクノイズ除去方法。 - Ｘ２≧Ｘ１−１である請求の範囲第１項乃至第４項のいずれか１項に記載のスパイクノイズ除去方法。
- 各測定点の位置は、受光素子のチャンネル位置である請求の範囲第１項乃至第５項のいずれか１項に記載のスパイクノイズ除去方法。
- 前記スパイクノイズ範囲の終点を判断するステップが、前記始点から所定の第１間隔Ｘ１の範囲内において、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが負値で｜Δｙｈ，ｈ＋１｜が、前記閾値（Ｙ１）以上にならないとき、当該範囲内に、スパイクノイズは存在しないとみなすステップ（Ｓ２５ａ）を有する請求の範囲第１項に記載のスパイクノイズ除去方法。
- 前記スパイクノイズ範囲の終点を判断するステップが、前記始点から所定の第１間隔Ｘ１の範囲内において、１次差分値又は１次微分値Δｙｈ，ｈ＋１＝ｙｈ＋１−ｙｈが正値でΔｙｈ，ｈ＋１が、前記閾値（Ｙ１）以上にならないとき、当該範囲内に、スパイクノイズは存在しないとみなすステップ（Ｓ２５ａ）を有する請求の範囲第２項に記載のスパイクノイズ除去方法。
- 請求の範囲第１項乃至第８項のいずれか１項に記載の計算方法をコンピュータに実行させるコンピュータプログラム。
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2003372677 | 2003-10-31 | ||
| JP2003372677 | 2003-10-31 | ||
| PCT/JP2004/010881 WO2005043085A1 (ja) | 2003-10-31 | 2004-07-23 | 平均化反復法を用いたスパイクノイズ除去方法及びコンピュータプログラム |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPWO2005043085A1 JPWO2005043085A1 (ja) | 2007-05-10 |
| JP4566913B2 true JP4566913B2 (ja) | 2010-10-20 |
Family
ID=34544044
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2005515089A Expired - Fee Related JP4566913B2 (ja) | 2003-10-31 | 2004-07-23 | 平均化反復法を用いたスパイクノイズ除去方法及びコンピュータプログラム |
Country Status (2)
| Country | Link |
|---|---|
| JP (1) | JP4566913B2 (ja) |
| WO (1) | WO2005043085A1 (ja) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20110301427A1 (en) * | 2010-06-04 | 2011-12-08 | Yongji Fu | Acoustic physiological monitoring device and large noise handling method for use thereon |
| JP6885995B2 (ja) * | 2019-09-13 | 2021-06-16 | 株式会社構造計画研究所 | 情報処理システム、情報処理方法及びプログラム |
| CN111259311B (zh) * | 2020-01-14 | 2023-03-24 | 西安应用光学研究所 | 尖峰噪声处理方法 |
| CN111274904B (zh) * | 2020-01-14 | 2023-03-24 | 西安应用光学研究所 | 针对尖峰噪声的信号处理系统 |
| CN111289489B (zh) * | 2020-03-05 | 2023-06-02 | 长春长光辰英生物科学仪器有限公司 | 一种基于拉曼光谱的微生物单细胞生长检测方法 |
| CN113376698B (zh) * | 2021-06-23 | 2022-09-27 | 吉林大学 | 一种低漏检的地震数据尖峰噪声检测及压制方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH06137889A (ja) * | 1992-10-22 | 1994-05-20 | Yamatake Honeywell Co Ltd | スパイクノイズフィルタ |
| JPH07103788A (ja) * | 1993-09-30 | 1995-04-18 | Nippon Seiko Kk | 信号弁別装置 |
| JPH08101042A (ja) * | 1994-09-29 | 1996-04-16 | Horiba Ltd | 分析計出力の処理方法 |
| JP2629511B2 (ja) * | 1991-12-20 | 1997-07-09 | 住友金属工業株式会社 | 信号中のノイズ除去方法 |
| JP2004020427A (ja) * | 2002-06-18 | 2004-01-22 | Saginomiya Seisakusho Inc | ノイズ除去方法及びノイズ除去フィルタ |
-
2004
- 2004-07-23 WO PCT/JP2004/010881 patent/WO2005043085A1/ja active Application Filing
- 2004-07-23 JP JP2005515089A patent/JP4566913B2/ja not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2629511B2 (ja) * | 1991-12-20 | 1997-07-09 | 住友金属工業株式会社 | 信号中のノイズ除去方法 |
| JPH06137889A (ja) * | 1992-10-22 | 1994-05-20 | Yamatake Honeywell Co Ltd | スパイクノイズフィルタ |
| JPH07103788A (ja) * | 1993-09-30 | 1995-04-18 | Nippon Seiko Kk | 信号弁別装置 |
| JPH08101042A (ja) * | 1994-09-29 | 1996-04-16 | Horiba Ltd | 分析計出力の処理方法 |
| JP2004020427A (ja) * | 2002-06-18 | 2004-01-22 | Saginomiya Seisakusho Inc | ノイズ除去方法及びノイズ除去フィルタ |
Also Published As
| Publication number | Publication date |
|---|---|
| JPWO2005043085A1 (ja) | 2007-05-10 |
| WO2005043085A1 (ja) | 2005-05-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US9159119B2 (en) | Method and system for super-resolution signal reconstruction | |
| US7379192B2 (en) | Optical metrology of single features | |
| US9784621B2 (en) | Spectroscopic apparatus and methods | |
| US11143555B2 (en) | Methods and devices for reducing spectral noise and spectrometry systems employing such devices | |
| US9163988B2 (en) | Detection systems and methods using coherent anti-stokes Raman spectroscopy | |
| US8712184B1 (en) | Method and system for filtering noises in an image scanned by charged particles | |
| WO2014064447A1 (en) | Spectroscopic apparatus and methods | |
| KR102120524B1 (ko) | 웨이퍼들 상의 결함들에 대한 정보의 결정 | |
| JP6023177B2 (ja) | 多項式フィッティングによるスペクトルデータのバックグラウンド放射線推定 | |
| JP4566913B2 (ja) | 平均化反復法を用いたスパイクノイズ除去方法及びコンピュータプログラム | |
| KR102514136B1 (ko) | 반도체 웨이퍼의 리피터 결함 포착 | |
| JP2007170941A (ja) | 歪み測定装置、方法、プログラムおよび記録媒体 | |
| JP2019502097A (ja) | 標的空間スペクトル検出を使用するタグ読み取り | |
| CN113837032B (zh) | 一种nv色心光探测磁共振曲线的极度欠采样重构方法 | |
| JP2007529721A (ja) | スペクトルデータ拡張方法 | |
| JP2016080668A (ja) | 表面処理状況モニタリング装置及び表面処理状況モニタリング方法 | |
| CN116754654B (zh) | 一种基于图像分析模型的声音屏障检测方法及装置 | |
| US5942763A (en) | Apparatus and method for identifying an identification mark of a wafer | |
| JP6560909B2 (ja) | プラズマ処理方法およびプラズマ処理装置 | |
| CN108051425B (zh) | 一种拉曼光谱信噪比评估方法 | |
| EP2520914A1 (en) | Estimation of background radiation in spectral data by polynomial fitting | |
| US7379186B2 (en) | Chirp indicator of ultrashort optical pulse | |
| US7973937B1 (en) | Near field suppression with a multi-aperture imaging system | |
| JP2007187459A (ja) | エッジ検出方法およびエッジ検出装置 | |
| WO2020203592A1 (ja) | 光検出装置、光検出方法、光検出装置の設計方法、試料分類方法、及び、不良検出方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20070621 |
|
| 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: 20100713 |
|
| 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: 20100804 |
|
| 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: 20130813 Year of fee payment: 3 |
|
| S531 | Written request for registration of change of domicile |
Free format text: JAPANESE INTERMEDIATE CODE: R313531 |
|
| S533 | Written request for registration of change of name |
Free format text: JAPANESE INTERMEDIATE CODE: R313533 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130813 Year of fee payment: 3 |
|
| 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 |
|
| LAPS | Cancellation because of no payment of annual fees |