JP3773563B2

JP3773563B2 - 2D digital signal processor

Info

Publication number: JP3773563B2
Application number: JP22289195A
Authority: JP
Inventors: 伸青木
Original assignee: Ricoh Co Ltd
Current assignee: Ricoh Co Ltd
Priority date: 1995-08-31
Filing date: 1995-08-31
Publication date: 2006-05-10
Anticipated expiration: 2015-08-31
Also published as: JPH0969755A

Description

【０００１】
【発明の属する技術分野】
本発明は、標本化によって生じる折り返し成分を打ち消し、原信号が持つ高周波成分を復元し、また、２次元デジタル信号を高解像度化する２次元デジタル信号処理装置に関する。
【０００２】
【従来の技術】
本出願人は先に、一次元のデジタル信号に対して、標本化位置を変えた複数の信号を使い、折り返し成分をキャンセルすることにより、折り返し歪を低減し、原信号が持つ高周波成分を復元する方法および装置を提案した（特願平７−１４２７７５号）。
【０００３】
【発明が解決しようとする課題】
上記した方法は、対象が一次元のデジタル信号か、あるいは標本化位置が縦（または横）方向だけにずれた２次元のデジタル信号に限定され、また、原信号がナイキスト周波数の２倍に帯域制限されていた。
【０００４】
しかし、カメラやスキャナなどデジタル画像入力機器では、その標本化位置のずれは一般に２次元的である。また、サンプリング前の原信号がナイキスト周波数の２倍の帯域に制限されているとは限らず、それ以上の高周波成分が含まれ、多重の折り返し歪が発生している場合もある。
【０００５】
本発明の目的は、２次元的な標本化位置のずれ、任意の数の折り返し成分を持つ複数の２次元デジタルデータに対して、その折り返し成分を打ち消すことにより、折り返し歪を低減し、原信号の高周波成分を復元する２次元デジタル信号処理装置を提供することにある。
【０００６】
【課題を解決するための手段】
前記目的を達成するために、請求項１記載の発明では、同一の信号について、標本化位置を変えて同一の標本化間隔でｎ回の標本化によって得られるｎ組の２次元デジタルデータに対して、標本化周波数の１／２以上の周波数帯域を含む、必要な原信号の帯域を全て透過するローパスフィルタ処理を行い、該ローパスフィルタ処理によって得られたｎ組の連続信号をそれぞれ前記標本化間隔より狭い間隔で標本化した画素値を高解像度化データとして求める高解像度化処理手段と、前記各デジタルデータの標本化位置に応じて、前記透過帯域に含まれる不要な折り返し成分を打ち消すような重みを算出する手段と、対応する位置における各高解像度化データに該重みをつけて加重和をとり出力する手段とを備えたことを特徴としている。
【０００７】
【実施例】
以下、本発明の一実施例を図面を用いて具体的に説明する。
まず、本発明で使用する信号処理法の原理について説明する。
標本化の対象となる２次元の原信号をｆ（ｘ，ｙ）とする。また、それをｘ，ｙについてそれぞれ間隔１で正方格子状に標本化した標本化信号をｄ（ｘ，ｙ）とする。
【０００８】
フーリエ変換、逆変換および畳み込みを式（１）、（２）、（３）のように定義する。
【０００９】
【数１】

【００１０】
〈標本化位置の異なるデジタルデータの周波数空間での形〉
標本化は、原信号と、以下のようなδパルス列（ｓ（ｘ，ｙ））との積をとる操作と考えられる。
【００１１】
【数２】

【００１２】
ここで、Ｓ（ｕ，ｖ）はｓ（ｘ，ｙ）のフーリエ変換で、周波数空間での表現である。以下、大文字の英字が周波数空間での表現であり、小文字の英字が実空間での表現である。
【００１３】
そうすると、標本化した信号（ｄ（ｘ，ｙ））およびそのフーリエ変換Ｄ（ｕ，ｖ）は、式（６）、（７）で表される。
【００１４】
【数３】

【００１５】
図２は、実空間での標本化の過程を説明する図であり、（ａ）はｆ（ｘ，ｙ）、（ｂ）はｓ（ｘ，ｙ）、（ｃ）はｆ（ｘ，ｙ）ｓ（ｘ，ｙ）である。ここで、ｆ（ｘ，ｙ）が画像信号であるとき、ｓ（ｘ，ｙ）が主、副走査信号となって、画像信号がサンプルされ（ｓ（ｘ，ｙ）の周期が画素間隔となる）、その振幅値（ｆ（ｘ，ｙ）ｓ（ｘ，ｙ））が量子化される。
【００１６】
図３は、実空間をフーリエ変換した周波数空間における標本化の過程を説明する図であり、（ａ）はｆ（ｘ，ｙ）のスペクトル、（ｂ）はｓ（ｘ，ｙ）のスペクトル、（ｃ）はｄ（ｘ，ｙ）のスペクトルであり、原信号ｆ（ｘ，ｙ）のフーリエ変換Ｆ（ｕ，ｖ）を間隔１の格子上に繰り返し並べた形になっている（標本化周波数＝１、ナイキスト周波数＝１／２）。以下、式（７）の（ｋ，ｌ）＝（０，０）の項を基本波、（ｋ，ｌ）＝（ｍ，ｎ）の項を（ｍ，ｎ）高調波と呼ぶ。
【００１７】
次に、同じｆ（ｘ，ｙ）について、標本化位置を（ｐ，ｑ）だけずらして同一の標本化間隔で標本化する場合を考える。この標本化した信号をｄｐ，ｑ（ｘ，ｙ）とする。これは、式（４）のδパルス列ｓ（ｘ，ｙ）を、（ｐ，ｑ）だけずらしたｓ（ｘ−ｐ，ｙ−ｑ）に変えたものである。位置をずらしたパルス列のフーリエ変換は、
【００１８】
【数４】

【００１９】
であり、これは、式（５）にｅｘｐ（−２πｉ（ｐｕ＋ｑｖ））を掛けた形となる。つまり、これはＳ（ｕ，ｖ）をｐｕ＋ｑｖに比例して位相を変化させた形になっており、基本波は式（５）のものと同じで、高調波はその位置に応じて位相がずれている。
【００２０】
そうすると、このずれた位置（ｐ，ｑ）で標本化した信号（ｄｐ，ｑ）のフーリエ変換は、
【００２１】
【数５】

【００２２】
となり、やはり基本波（（ｋ，ｌ）＝（０，０））は式（７）と同じ（Ｆ（ｕ，ｖ））で、高調波はその中心位置に応じた分だけ位相を変えた形になる。
【００２３】
このように、折り返し歪を含み、標本化位置の異なる複数の標本化信号を周波数空間で見ると、基本波は標本化位置（ｐ，ｑ）によらず原信号と同一で、（ｋ，ｌ）高調波は、標本化位置に応じて−２π（ｐｋ＋ｑｌ）だけ位相がずれていることが分かる。この位相のずれを使って、以下説明するように、折り返された成分だけを取り除き、本来の高周波成分を復元することができる。
【００２４】
〈高周波成分の復元〉
以上の解析結果を使い、一つの原信号を異なる位置で標本化した複数の組の標本化信号から、折り返し成分を取り除き、原信号の高周波成分を復元する原理を説明する。
【００２５】
ここでは、原信号ｆ（ｘ，ｙ）が、周波数空間上で、ナイキスト周波数（１／２）の２倍に帯域制限されているものとする。つまり、Ｆ（ｕ，ｖ）は、図４に示すように、ｕ＜−１，１＜ｕ，ｖ＜−１，１＜ｖの範囲では、その値が０となるものとする。
【００２６】
原信号がこのように帯域制限されていれば、標本化信号は、どの周波数成分も基本波、高調波の内、たかだか２つを足したものとなる（図５）。
【００２７】
また、ここでは標本化信号の標本化位置は既知であると仮定し、位置（ｐ０，ｑ０），（ｐ１，ｑ１）...で標本化した標本化信号をそれぞれｄ０（ｘ，ｙ）、ｄ１（ｘ，ｙ）...とし、またそれぞれのフーリエ変換をＤ０（ｕ，ｖ）、Ｄ１（ｕ，ｖ）...とする。
【００２８】
まず、一つの標本化信号ｄ０（ｘ，ｙ）にローパスフィルタをかけ、（−１＜ｕ，ｖ＜１）の帯域だけを取り出す。この信号をｅ０（ｘ，ｙ）とする。そのフーリエ変換（Ｅ０（ｕ，ｖ））は次のようになる。
【００２９】
【数６】

【００３０】
また、Σ_k，_lは、−１，０，１について、つまり、
（ｋ，ｌ）＝（−１，−１），（−１，０），（−１，１），（０，−１），（０，０），（０，１），（１，−１），（１，０），（１，１）の９通りについての和をとる。これは、図６に示すように、基本波と周波数空間で上下左右斜めに隣接する高調波の一部が重なった信号である。
【００３１】
同様に、その他の標本化信号ｄ１（ｘ，ｙ）...もそれぞれＬＰＦ処理し、ｅ１（ｘ，ｙ）、ｅ２（ｘ，ｙ）...を求める。これら信号は、式（９）より、基本波は同一で、高調波は標本化位置（（ｐ０，ｑ０），（ｐ１，ｑ１）...）に応じて、その位相がずれている。
【００３２】
次に、これら複数の信号を、適当な重みをつけて加重和をとる。ｗｎで加重和をとった信号、
【００３３】
【数７】

【００３４】
は、周波数空間で表現すれば、式（１１）より、
【００３５】
【数８】

【００３６】
となる。
【００３７】
よって、式（１４）から式（２１）までの８つの項が０となり、式（１３）の項だけが残れば、（−１＜ｕ，ｖ＜１）の範囲の原信号が完全に復元できることになる。
【００３８】
そのためには重み（ｗｎ）が以下の条件を満たせばよい。まず、Ｆ（ｕ，ｖ）をそのまま残すため式（１３）より
【００３９】
【数９】

【００４０】
なる条件が必要である。
【００４１】
また、式（１４）の項が０となるためには、
【００４２】
【数１０】

【００４３】
なる条件を満たせばよい。対称性により、この条件で式（２１）の項も０となる。その他の項も同様であり、式（２３）も含めてまとめると、（ｋ，ｌ）＝（−１，−１），（−１，０），（−１，１），（０，−１）の４通りについて、
【００４４】
【数１１】

【００４５】
を満たせば、式（１４）から式（２１）の８つの項が全て０となる。
【００４６】
以上の式（２２）および式（２６）の条件を満たす重み（ｗｎ）を使用すれば、式（１２）は結局
【００４７】
【数１２】

【００４８】
となり、原信号（Ｆ（ｕ，ｖ)）を復元することができる。式（２２）、(２７）、（２８）は、重みｗｎについての連立１次方程式であり、式（２２）が１本、式（２７)）（２８）が４通りの（ｋ，ｌ）について８本、計９本の連立方程式なので、９つの異なる位置で標本化したデジタルデータ（ｄ０（ｘ，ｙ）...ｄ８（ｘ，ｙ））があれば、条件を満たす解（ｗｎ）を求めることができる。また、これは連立１次方程式なので、行列演算によって簡単に解くことができる。
【００４９】
〈様々な帯域の原信号への対応〉
以上は、原信号（ｆ（ｘ，ｙ））が、ナイキスト周波数の２倍（−１＜ｕ，ｖ＜１）に帯域制限されているものと仮定して説明した。原信号がそれ以外の帯域を持つ場合でも、同様に折り返し成分を打ち消すことができる。
【００５０】
例えば、原信号が横方向にはナイキスト周波数の３倍に、縦方向には２倍に制限されている場合、基本波に重なる高調波は、図７に示すようになる。この場合では、原信号の帯域に重なる高調波は、（ｋ＝−２，−１，０，１，２），（ｌ＝−１，０，１）の組み合わせで、計１４個ある。
【００５１】
そこで、式（１０）のＬＰＦの代わりに、この原信号を全て透過するようなＬＰＦつまり、
【００５２】
【数１３】

【００５３】
を用いてフィルタ処理すると、得られる信号（ｅｎ）は、式（１２）の（１４）から（２１）の８つの項に、（ｋ＝−２，２）に対応する６つの項が加わり、高調波成分が合計で１４項のものになる。
【００５４】
それらを全て打ち消すためには、対称性を考慮して７通りの（ｋ，ｌ）の組み合わせについて、１４個の方程式からなる式（２７）、（２８）の条件を満たせばよい。この場合、重みｗｎを求めるには、式（２２）と合わせて１５本の連立方程式を解くため、１５組の入力データが必要となる。
【００５５】
このように、原信号の復元すべき帯域に応じて、必要な帯域を全て透過するようなＬＰＦの帯域を選択し、その帯域に重なる不要な高調波を全て打ち消すように、重みｗｎを決定することによって、原信号を復元することができる。一般にＭ個の折り返し成分を打ち消すためには、２Ｍ＋１組の入力データが必要となる。
【００５６】
〈雑音への対応〉
上記した説明では、各信号には雑音は含まれないものと考え、原信号の完全な復元を行うものであった。しかし、現実には信号に必ず雑音が含まれるので、雑音に強い信号処理方法が求められる。
【００５７】
ここでは、まず前述したような加重和による信号の雑音の現れ方を解析し、それに対する２つの対処方法を説明する。
【００５８】
一般に、複数の信号が、互いに無相関で分散σ²の雑音を含むとすると、重みｗｎを使ったそれらの加重平均信号の雑音は、
【００５９】
【数１４】

【００６０】
なる分散を持つ。つまり、加重平均信号の雑音は、それぞれの信号のΣ（ｗｎ）²倍になると考えられる。
【００６１】
式（２２）、（２６）で求められる重み（ｗｎ）は、雑音を考慮していないので、標本化位置の組み合わせによっては、Σ（ｗｎ）²が非常に大きくなる場合があり、得られる加重平均信号は折り返し成分が打ち消されると同時に雑音が強調される可能性がある。そこで、雑音の低減も考慮した重みを考える。
【００６２】
〈「折り返し成分を打ち消し」かつ「雑音最小化」のための重み〉
９つの入力信号を使う場合は、折り返し成分を打ち消すような加重（ｗｎ）は、式（２６）の条件から一意に決まる。一方、１０組以上の入力信号を使えば、自由度が増えるので、式（２６）を満たした上で、残る自由度を使って式(３０)を最小化する重みをとることにより、折り返し成分を打ち消した上で、更に雑音の低減を図ることができる。
【００６３】
つまり、式（２２）、（２６）を満たす重みの代わりに、
【００６４】
【数１５】

【００６５】
なる条件のもとで、
【００６６】
【数１６】

【００６７】
を最小化する重みを使う。
【００６８】
これはｗｎについて、１次式の条件付き、２乗和の最小化であるので、公知の計算法（Ｌａｇｒａｎｇｅの未定乗数法）により、連立１次方程式を解くことで簡単に解が得られる。
【００６９】
〈「折り返し成分および雑音」の最小化のための重み〉
雑音成分との兼ね合いを考えると、折り返し成分も一種の雑音であり、必ずしも０にする必要はない。折り返し成分および雑音成分をそれぞれ低減する重みも有効である。
【００７０】
これは、式（２６）を束縛条件とはせず、次式のような、その絶対値の２乗値
【００７１】
【数１７】

【００７２】
を式（３０）と同等の雑音評価値と考え、それぞれの雑音評価値に対して増加関数となる評価関数をとり、それを最小化する重みを採用することによって実現できる。
【００７３】
Σｗｎ＝１の条件のもとで、例えば次式のような、式（３２）と式（３０）の加重和を評価関数とし、
【００７４】
【数１８】

【００７５】
これを最小化する重みを使う。この場合、α≒０では折り返し成分の打ち消しが優先となり、α≒１では基本波の雑音低減が優先となる。式（３３）も、ｗｎについての２乗和の最小化であり、やはり、連立１次方程式を解くことで簡単に解が得られる。この場合は、２組以上の入力データが必要である。
【００７６】
〈標本化位置の推定〉
上記した説明では、それぞれの入力標本化信号について、標本化位置は既知であるとして考えてきた。しかし、デジタルカメラを動かしながら撮影したデータを扱う場合などにおいては、標本化位置を予め知ることができないので、与えられたデジタルデータからそれぞれの標本化位置を推定する必要がある。その推定方法を以下、説明する（なお、ここで用いるアルゴリズムは、「画像の時空間微分算法を用いた速度ベクトルの分布計測システム」（計測自動制御学会論文集）Ｖｏｌ．２２，Ｎｏ．１２，ｐｐ８８−９４（昭和６１年１２月）に記載されたものである）。
【００７７】
２つのデジタルデータの標本化位置のずれは、移動する対象物の速度とも考えることができる。そこで、速度（ｓ，ｔ）で平行移動する１次元の連続関数ｆ（ｘ，ｙ，ｔ）を例に考える。変化が微小ならば、Ｔａｙｌｅｒ展開の１次近似として、
【００７８】
【数１９】

【００７９】
が各点で成り立つ。この式の∂ｆ／∂ｘ，∂ｆ／∂ｙおよび∂ｆ／∂ｔから、各点で変化速度（ｓ，ｔ）を求めることができる。
【００８０】
しかし、上記した式では、
（１）１次近似からの外れ
（２）雑音の影響
などのため精度が悪いので、次式のような、式（３４）の左辺の２乗をある領域内で積分した値を評価関数とし、この評価関数が最小値をとる（ｓ，ｔ）を求めることで、正確な値が得られる。
【００８１】
【数２０】

【００８２】
このようにして、２つのデータ間の相対位置が分かれば、１つのデータを基準として、他の全てのデータとの相対位置を求めることにより、全てのデータの相対位置を推定することができる。
【００８３】
〈実施例〉
図１は、本発明の実施例のブロック構成図である。本実施例では、前述した本発明の原理に基づき、同一の対象を異なる複数の標本化位置で標本化した２次元デジタルデータを入力し、標本化される前の原信号の高周波成分を復元した、高解像度のデジタルデータを得るための信号処理方法を説明する。以下、各部分について、その機能を説明する。
【００８４】
（入力デジタルデータ）
本発明で使用する入力デジタルデータ１は、同一の対象を、位置をずらして同一の標本化間隔で標本化した２次元デジタルデータである。例えば、デジタルカメラを使い、同一の対象物を少しずつ移動しながら撮影した複数のデータや、ビデオカメラで同一の対象物をゆっくりと移動しながら撮影したシーンの連続する複数のコマをデジタル化したデータなどが利用できる。また、カメラなど入力装置の光学系のＭＴＦやセンサの開口特性を予め調べることにより、標本化される前の原信号の帯域が分かる。ここでは、ナイキスト周波数の２倍に帯域制限された原信号を復元するための例を示す。
【００８５】
以下、入力デジタルデータをそれぞれｉｎ０，ｉｎ１，...と呼ぶ。また、第ｎ組の入力デジタルデータ（ｉｎ_n）の、第ｉ行、第ｊ列のサンプル（画素）を、ｉｎ_n〔ｉ，ｊ〕と呼ぶ。
【００８６】
（基準データバッファ）
基準データバッファ２は、標本化位置の基準の位置とするため、第０組のデータ（ｉｎ０）を格納するバッファである。後段の位置差推定部では、２組のデジタルデータを入力し、それらの標本化位置の差を推定する。第０組のデータ（ｉｎ０）を基準データとし、ｉｎ０とｉｎ１、およびｉｎ０とｉｎ２の標本化位置差を推定することにより、第０組のデータの標本化位置を基準として、全ての入力データの相対的な標本化位置を推定することができる。
【００８７】
（位置差推定部）
位置差推定部３は、それぞれの入力データ（ｉｎ１，ｉｎ２...）について、基準データバッファに格納された基準データ（ｉｎ０）との、標本化位置の差を推定することにより、各入力データの相対的な標本化位置を推定する手段である。
【００８８】
第１組のデータの標本化位置を（ｐ１，ｑ１）、第２組のデータの標本化位置を（ｐ２，ｑ２）などとして出力する（基準となる第０組のデータの標本化位置は（０，０）である）。このような推定は、式（３５）の微分、積分をデジタルデータの差分、加算に置き換えることにより実現できる。
【００８９】
すなわち、ｉｎ０とｉｎ_nの位置差（ｐｉ，ｑｉ）を推定するには、式（３５）を
【００９０】
【数２１】

【００９１】
置き換え、また、対象の速度（ｓ，ｔ）ではなく、標本化位置の移動量を測るので、符号を変えて、
【００９２】
【数２２】

【００９３】
を最小化する（ｐｎ，ｑｎ）として求められる。これは線形最小２乗法により解くことができる。
【００９４】
（加重計算部）
加重計算部４は、位置差推定部３で得られた標本化位置（（０，０），（ｐ１，ｑ１）...）に基づいて、後段の積和部で使用される加重を求める手段である。第０，１...組の入力デジタルデータに対する加重として、それぞれｗ０，ｗ１，...を出力する。
【００９５】
具体的には、ここでは原信号はナイキスト周波数の２倍の帯域を持つので、原理の説明で述べたように９組の入力デジタルデータを使い、（ｋ，ｌ＝−１，０，１）として、式（２２）、（２６）を解く演算を行う。また、他の方法として、１０組以上の入力デジタルデータを使い、式（２２）、（２６）の束縛条件のもとで、式（３１）を最小化する重みを求めることもできる。さらに、他の方法として、２組以上の入力デジタルデータを使い、式（２２）の束縛条件のもとで、式（３３）を最小化する重みを求めることもできる。また、原信号の帯域がこれと異なる場合には、不要な（ｋ，ｌ）高調波を全て打ち消すような重みｗｎを求める。
【００９６】
（解像度倍率指定部）
本実施例の最終的な出力デジタルデータは、入力デジタルデータよりも画素数の多い高解像度なデータである。解像度倍率指定部５は、出力デジタルデータの解像度を、入力デジタルデータの解像度の倍率として指定する手段である。ここでは、解像度を４倍にするものとする。つまり、出力デジタルデータは入力デジタルデータに比べ、サンプル数が４倍、標本化周期が１／４となる。
【００９７】
解像度の倍率は、利用者が指定するようにしてもよい。ただし、本実施例では、ナイキスト周波数の２倍の周波数成分までを復元するので、出力データに折り返し歪を含めないためには、倍率は２倍以上としなければならない。
【００９８】
（広帯域ＬＰＦ処理および高解像度化処理部）
この部分は、本発明の原理で説明した「〔−１，１〕の透過帯域を持つＬＰＦにより、基本波全部および隣接する高調波の一部を足した信号を取り出す」部分に相当する。
【００９９】
入力デジタルデータは、ｉｎ_n〔ｉ，ｊ〕というデータ列であるが、これを次式のように、位置（ｘ，ｙ）に対する連続関数（デルタパルス列）と考えたものが式（６）の標本化信号ｄ（ｘ，ｙ）にあたる。
【０１００】
【数２３】

【０１０１】
この信号（ｄｎ（ｘ，ｙ））に、式（１０）の広帯域ＬＰＦを掛ける。
【０１０２】
理想的なＬＰＦ処理は、周波数空間では、
【０１０３】
【数２４】

【０１０４】
なる矩形関数との乗算であり、実空間では、
【０１０５】
【数２５】

【０１０６】
とのコンボリューションである。ｄｎ（ｘ，ｙ）をＬＰＦ処理した信号をｅｎ（ｘ，ｙ）とすると、
【０１０７】
【数２６】

【０１０８】
である。
【０１０９】
さらに、式（３７）より、入力デジタルデータ（ｉｎ_n〔ｉ，ｊ〕）を使って書き換えると、
【０１１０】
【数２７】

【０１１１】
と書ける。
【０１１２】
このＬＰＦ処理によって得られる信号（ｅｎ（ｘ，ｙ））は連続信号であり、これを標本化した信号を、このブロック６の出力デジタル信号とする。そのため、連続信号の内、欲しい（標本化される）位置（ｘ，ｙ）についてだけ式(３８)の計算をすればよい。
【０１１３】
ここで得られるデジタル信号は、後段の積和部８で加算される。加算は、原信号に対して同じ位置同志の画素値を足しあわせるものである。そのためｅｎを標本化する位置は、それぞれ標本化位置の異なる各入力画像上での、原画像上の同じ位置に対応する位置とするのがよい。
【０１１４】
ここでは、基準データｉｎ０〔ｉ，ｊ〕の高解像度化データｍｉｄ０〔ｉ，ｊ〕の画素位置を基準とし、その他の各入力データについては、高解像度化基準データの画素位置に対応する、各画像上の位置の画素値を求める。
【０１１５】
図８は、高解像度化画素の位置関係を示す。図中、大きな白丸○は、標本化間隔１で標本化された基準入力データの画素を示す。大きな黒丸●は、基準入力データに対して（ｐｎ，ｑｎ）ずれた位置で、標本化間隔１で標本化された第ｎ組の入力データの画素位置を示す。
【０１１６】
高解像度化処理部は、まず、基準入力データに対しては、単純に標本化間隔を１／４とし、位置（ｉ／４，ｊ／４）（ｉ，ｊは整数）の画素値を求める。図８において、小さな白丸が基準入力データを高解像度化した画素位置を示す。
【０１１７】
第ｎ組の入力データ（ｉｎ_n）の高解像度化データは、次のように作成する。
すなわち、第ｎ組の入力データの画素値（黒丸）を使用して（例えば、補間する）、原信号上で高解像度化基準データの画素位置（図８の小さな白丸）と同じ位置にあたる画素値を第ｎ組の高解像度化データとして求める。図９は、第ｎ組の入力データの画素値を例えば補間することによって、作成された第ｎ組の高解像度化データの画素位置を示す。図の小さな黒丸の位置は、基準入力データを高解像度化した画素位置（小さな白丸）と同じ位置になる。第ｎ組の入力データは、標本化位置が（ｐｎ，ｑｎ）だけずれているので、（（ｉ／４）−ｐｎ，（ｊ／４）−ｑｎ）が対応する位置となる。
【０１１８】
なお、先に提案した方法では、高解像度化データの画素位置は、各入力データに相対的に決め、それらの加重和をとるときに位置をずらして加算した。この方法では、加算するために画素位置をずらす単位は、高解像度化データの画素間隔となる。例えば、高解像度化データの倍率を４倍とすると、入力データの１／４画素単位でしかずらすことができない。そのため、倍率が小さい場合には、位置合わせの精度が悪く、結果として、復元される信号の精度が悪くなる恐れがあった。これに対して、本発明の方法は、第ｎ組の入力データ（つまり画素値）と、その標本化位置（ｐｎ，ｑｎ）を基に、基準となる高解像度化データの画素位置と同じ位置の画素値を求めているので、高解像度化の倍率によらず、高解像度化データの画素位置を高精度に決定することができる。
【０１１９】
このようにして得られる「広帯域ＬＰＦ処理と高解像度化処理部」の出力となるデジタルデータをｍｉｄｎ〔ｉ，ｊ〕と呼ぶ。ここでは新たな標本化間隔は１／４としたので、ｍｉｄｎ〔ｉ，ｊ〕は結局、
【０１２０】
【数２８】

【０１２１】
となり、予めｓｉｎｃ２（ｉ，ｊ）なるテーブル（実際には有限の範囲内でそれを近似したもの）を用意しておけば、入力デジタル信号ｉｎ_nの積和によって計算できる。
【０１２２】
ＬＰＦ処理された信号ｅｎ（ｘ，ｙ）は、原信号のナイキスト周波数よりも高い周波数成分を含む信号である。そしてそれを有効に表現するためには、つまり原信号の高周波成分を復元するためには、ＬＰＦ処理された信号ｅｎ（ｘ，ｙ）に対するサンプリング間隔は、入力信号のサンプリング間隔より狭い間隔でなければならない。信号ｅｎ（ｘ，ｙ）に対するサンプリング間隔が、入力信号のサンプリング間隔より狭くなければ、再び折り返し歪が発生し、意味がない。そのため、ここでは１／４としている。また、原信号の帯域がこれと異なる場合には、原信号を復元するのに必要な帯域を透過するようなＬＰＦを用いる。
【０１２３】
（中間データバッファ）
中間データバッファ７は、前段の「広帯域ＬＰＦ処理と高解像度化処理部」で得られた中間デジタルデータ（ｍｉｄｎ〔ｉ，ｊ〕）を保存するバッファである。
【０１２４】
（積和部）
この積和部７での加重和計算は、本発明の原理で説明した「式（１２）の加重和をとり、高調波成分を打ち消す」部分に相当する。本実施例の構成では、実空間で処理する。実空間でもやはり各中間データに、それぞれ加重ｗｎをかけた後、加算すればよい。
【０１２５】
中間データバッファ７には、原信号上で同じ位置に相当するデータが格納されているので、それぞれのバッファの画素値（ｍｉｄ０，ｍｉｄ１...）に重みｗ０，ｗ１...を掛けて足しあわせればよい。
【０１２６】
最終的に出力されるデータ（ｏｕｔ〔ｉ，ｊ〕）は、
【０１２７】
【数２９】

【０１２８】
となる。
【０１２９】
【発明の効果】
以上、説明したように、本発明によれば、同一の信号を、標本化位置を変えて同一の標本化間隔でｎ回の標本化によって得られるｎ組の２次元デジタルデータについて、必要な原信号の帯域を全て透過するフィルタと、該フィルタを透過した不要な高調波を全て打ち消す重みを使った加重和をとることにより、標本化による不要な折り返し歪を全て打ち消し、必要な原信号の高周波成分をすべて復元することができる。これにより、ぼけ、歪の少ない高解像度のデータを得ることが可能となる。
【０１３０】
また、本発明によれば、雑音を含まないデータに対して、折り返し成分を完全に打ち消すので、高精度に原信号の高周波成分を復元することができる。
【０１３１】
また、本発明によれば、雑音を含むデータに対しても、折り返し成分と雑音を同時に低減することができる。
【０１３２】
また、本発明によれば、高解像度化の倍率によらず、正確に折り返し成分を打ち消し、原信号の高周波成分を復元することができる。
【図面の簡単な説明】
【図１】本発明の実施例のブロック構成図である。
【図２】（ａ）、（ｂ）、（ｃ）は、実空間での標本化の過程を説明する図である。
【図３】（ａ）、（ｂ）、（ｃ）は、周波数空間における標本化の過程を説明する図である。
【図４】ナイキスト周波数の２倍に帯域制限された原信号を示す。
【図５】標本化により折り返しが起こった信号を示す。
【図６】ＬＰＦ処理された信号を示す。
【図７】標本化により３重の折り返しが起こった信号を示す。
【図８】高解像度化画素の位置関係を示す。
【図９】第ｎ組の高解像度化画素の位置を示す。
【符号の説明】
１入力デジタルデータ
２基準データバッファ
３位置差推定部
４加重計算部
５出力解像度指定部
６広帯域ＬＰＦ処理および高解像度化処理部
７中間データバッファ
８積和部
９出力デジタルデータ[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a two-dimensional digital signal processing apparatus that cancels aliasing components generated by sampling, restores high-frequency components of an original signal, and increases the resolution of a two-dimensional digital signal.
[0002]
[Prior art]
The applicant first uses a plurality of signals with different sampling positions for a one-dimensional digital signal, cancels the aliasing component, reduces aliasing distortion, and restores the high-frequency component of the original signal. Proposed a method and apparatus (Japanese Patent Application No. 7-142775).
[0003]
[Problems to be solved by the invention]
The above-described method is limited to a one-dimensional digital signal or a two-dimensional digital signal whose sampling position is shifted only in the vertical (or horizontal) direction, and the original signal has a bandwidth of twice the Nyquist frequency. It was restricted.
[0004]
However, in a digital image input device such as a camera or a scanner, the sampling position deviation is generally two-dimensional. Further, the original signal before sampling is not necessarily limited to a band twice the Nyquist frequency, and there may be a case where multiple higher-frequency components are included and multiple aliasing distortion occurs.
[0005]
An object of the present invention is to reduce aliasing distortion by canceling aliasing components of a plurality of two-dimensional digital data having a two-dimensional sampling position deviation and an arbitrary number of aliasing components, thereby reducing the original signal. Another object of the present invention is to provide a two-dimensional digital signal processing apparatus that restores the high-frequency component of the signal.
[0006]
[Means for Solving the Problems]
To achieve the above object, according to the first aspect of the present invention, n sets of two-dimensional digital data obtained by sampling n times at the same sampling interval while changing the sampling position for the same signal. Low-pass that transmits all necessary original signal bands, including a frequency band of 1/2 or more of the sampling frequency. Perform filtering, The Low pass Filtering Obtained by n pairs Continuous Each signal Obtain pixel values sampled at intervals narrower than the sampling interval as high resolution data Higher resolution processing Depending on the means and the sampling position of each digital data, Said A means for calculating a weight for canceling an unnecessary aliasing component included in the transmission band, and a means for adding a weight to each resolution-enhanced data at a corresponding position and outputting the weighted sum are provided. Yes.
[0007]
【Example】
Hereinafter, an embodiment of the present invention will be described in detail with reference to the drawings.
First, the principle of the signal processing method used in the present invention will be described.
A two-dimensional original signal to be sampled is assumed to be f (x, y). Further, a sampling signal obtained by sampling it in a square lattice pattern at intervals of 1 for x and y is assumed to be d (x, y).
[0008]
Fourier transform, inverse transform, and convolution are defined as shown in equations (1), (2), and (3).
[0009]
[Expression 1]

[0010]
<Shape of digital data at different sampling positions in frequency space>
Sampling is considered to be an operation of taking the product of the original signal and the following δ pulse train (s (x, y)).
[0011]
[Expression 2]

[0012]
Here, S (u, v) is a Fourier transform of s (x, y) and is an expression in the frequency space. In the following, uppercase alphabetic characters are representations in the frequency space, and lowercase alphabetic characters are representations in the real space.
[0013]
Then, the sampled signal (d (x, y)) and its Fourier transform D (u, v) are expressed by equations (6) and (7).
[0014]
[Equation 3]

[0015]
FIG. 2 is a diagram for explaining the sampling process in real space, where (a) is f (x, y), (b) is s (x, y), and (c) is f (x, y). ) S (x, y). Here, when f (x, y) is an image signal, s (x, y) is a main and sub-scan signal, and the image signal is sampled (the period of s (x, y) is the pixel interval. The amplitude value (f (x, y) s (x, y)) is quantized.
[0016]
FIG. 3 is a diagram for explaining a sampling process in a frequency space obtained by Fourier transforming a real space, where (a) is a spectrum of f (x, y), (b) is a spectrum of s (x, y), (C) is a spectrum of d (x, y), and has a form in which Fourier transforms F (u, v) of the original signal f (x, y) are repeatedly arranged on a grid with an interval of 1 (sampling). Frequency = 1, Nyquist frequency = 1/2). In the following, the term (k, l) = (0, 0) in equation (7) is called the fundamental wave, and the term (k, l) = (m, n) is called the (m, n) harmonic.
[0017]
Next, consider a case where the same f (x, y) is sampled at the same sampling interval by shifting the sampling position by (p, q). Let this sampled signal be dp, q (x, y). This is obtained by changing the δ pulse train s (x, y) in Expression (4) to s (xp, yq) shifted by (p, q). The Fourier transform of the shifted pulse train is
[0018]
[Expression 4]

[0019]
This is a form obtained by multiplying expression (5) by exp (-2πi (pu + qv)). In other words, this is a form in which the phase is changed in proportion to pu + qv of S (u, v), the fundamental wave is the same as that of Equation (5), and the harmonic wave has a phase depending on its position. It's off.
[0020]
Then, the Fourier transform of the signal (dp, q) sampled at this shifted position (p, q) is
[0021]
[Equation 5]

[0022]
The fundamental wave ((k, l) = (0,0)) is the same as (7) (F (u, v)), and the phase of the harmonic wave is changed by the amount corresponding to the center position. Become a shape.
[0023]
As described above, when a plurality of sampling signals including aliasing distortion and having different sampling positions are viewed in the frequency space, the fundamental wave is the same as the original signal regardless of the sampling positions (p, q), and (k, l ) It can be seen that the harmonics are out of phase by -2π (pk + ql) depending on the sampling position. Using this phase shift, as described below, only the folded components can be removed and the original high frequency components can be restored.
[0024]
<Restoring high frequency components>
Using the above analysis results, the principle of removing the aliasing component from a plurality of sets of sampled signals obtained by sampling one original signal at different positions and restoring the high frequency component of the original signal will be described.
[0025]
Here, it is assumed that the original signal f (x, y) is band-limited to twice the Nyquist frequency (1/2) in the frequency space. That is, F (u, v) is assumed to have a value of 0 in the range of u <-1, 1 <u, v <-1, 1 <v, as shown in FIG.
[0026]
If the original signal is band-limited in this way, the sampled signal is obtained by adding at most two of the fundamental wave and the harmonic wave to any frequency component (FIG. 5).
[0027]
Here, assuming that the sampling position of the sampling signal is known, the sampling signals sampled at the positions (p0, q0), (p1, q1)... Are respectively d0 (x, y), d1 (x, y)... and the respective Fourier transforms are D0 (u, v), D1 (u, v).
[0028]
First, a low-pass filter is applied to one sampled signal d0 (x, y) to extract only the band of (−1 <u, v <1). This signal is assumed to be e0 (x, y). The Fourier transform (E0 (u, v)) is as follows.
[0029]
[Formula 6]

[0030]
Also, Σ _k , _l Is for -1, 0, 1, ie
(K, l) = (-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 0), (0, 1), (1, −1), (1, 0), and (1, 1) are summed for nine patterns. As shown in FIG. 6, this is a signal in which a part of harmonics adjacent to each other in the vertical and horizontal directions in the frequency space overlap with the fundamental wave.
[0031]
Similarly, the other sampling signals d1 (x, y)... Are also subjected to LPF processing to obtain e1 (x, y), e2 (x, y). These signals have the same fundamental wave and the phase of the harmonic wave is shifted according to the sampling position ((p0, q0), (p1, q1)...) According to equation (9).
[0032]
Next, a weighted sum is obtained by applying appropriate weights to the plurality of signals. a signal with a weighted sum of wn,
[0033]
[Expression 7]

[0034]
Is expressed in frequency space, from Equation (11),
[0035]
[Equation 8]

[0036]
It becomes.
[0037]
Therefore, if the eight terms from Equation (14) to Equation (21) are 0 and only the term of Equation (13) remains, the original signal in the range of (−1 <u, v <1) is completely restored. It will be possible.
[0038]
For that purpose, the weight (wn) should satisfy the following conditions. First, from equation (13) to leave F (u, v) as it is
[0039]
[Equation 9]

[0040]
The following conditions are necessary.
[0041]
In order for the term of equation (14) to be 0,
[0042]
[Expression 10]

[0043]
It is sufficient to satisfy the following conditions. Due to symmetry, the term in equation (21) is also zero under this condition. The other terms are also the same. When the formula (23) is also included, (k, l) = (− 1, −1), (−1, 0), (−1, 1), (0, − About 4 ways of 1)
[0044]
[Expression 11]

[0045]
If all are satisfied, all the eight terms in the equations (14) to (21) become zero.
[0046]
If the weight (wn) satisfying the conditions of the above equations (22) and (26) is used, the equation (12) is eventually
[0047]
[Expression 12]

[0048]
Thus, the original signal (F (u, v)) can be restored. Equations (22), (27), and (28) are simultaneous linear equations with respect to the weight wn, with one equation (22) and four (k, l) equations (27) and (28). Since there are a total of 9 simultaneous equations, there are digital data (d0 (x, y)... D8 (x, y)) sampled at 9 different positions, so that the solution satisfying the condition (wn) Can be requested. Further, since this is a simultaneous linear equation, it can be easily solved by matrix operation.
[0049]
<Support for original signals in various bands>
The above has been described on the assumption that the original signal (f (x, y)) is band-limited to twice the Nyquist frequency (−1 <u, v <1). Even when the original signal has a band other than that, the aliasing component can be canceled in the same manner.
[0050]
For example, when the original signal is limited to three times the Nyquist frequency in the horizontal direction and twice in the vertical direction, the harmonics overlapping the fundamental wave are as shown in FIG. In this case, there are a total of 14 harmonics in a combination of (k = −2, −1, 0, 1, 2) and (l = −1, 0, 1) overlapping the band of the original signal.
[0051]
Therefore, instead of the LPF of Equation (10), an LPF that transmits all of the original signal, that is,
[0052]
[Formula 13]

[0053]
When the signal (en) is filtered using, the six terms corresponding to (k = −2, 2) are added to the eight terms (14) to (21) of the equation (12), The total number of harmonic components is 14 items.
[0054]
In order to cancel all of them, it is only necessary to satisfy the conditions of the equations (27) and (28) consisting of 14 equations for seven (k, l) combinations in consideration of symmetry. In this case, in order to obtain the weight wn, 15 simultaneous equations are solved together with the equation (22), so 15 sets of input data are required.
[0055]
In this way, according to the band to be restored of the original signal, an LPF band that transmits all necessary bands is selected, and the weight wn is determined so as to cancel all unnecessary harmonics that overlap the band. Thus, the original signal can be restored. In general, 2M + 1 sets of input data are required to cancel M folding components.
[0056]
<Response to noise>
In the above description, it is assumed that each signal does not contain noise, and the original signal is completely restored. However, in reality, since a signal always includes noise, a signal processing method that is resistant to noise is required.
[0057]
Here, first, the appearance of signal noise due to the weighted sum as described above is analyzed, and two coping methods are explained.
[0058]
In general, multiple signals are uncorrelated with each other and distributed σ ² , The noise of those weighted average signals using the weight wn is
[0059]
[Expression 14]

[0060]
With a variance of That is, the noise of the weighted average signal is Σ (wn) of each signal. ² It will be doubled.
[0061]
Since the weight (wn) obtained by the equations (22) and (26) does not consider noise, Σ (wn) depends on the combination of sampling positions. ² May become very large, and the resulting weighted average signal may cancel out the aliasing component and simultaneously enhance noise. Therefore, a weight considering noise reduction is considered.
[0062]
<Weights for "Canceling aliasing" and "Minimizing noise">
When nine input signals are used, the weight (wn) for canceling the aliasing component is uniquely determined from the condition of Expression (26). On the other hand, if more than 10 sets of input signals are used, the degree of freedom increases. Therefore, by satisfying Equation (26) and taking the weight to minimize Equation (30) using the remaining degree of freedom, the aliasing component is obtained. In addition, the noise can be further reduced.
[0063]
That is, instead of the weights that satisfy Equations (22) and (26),
[0064]
[Expression 15]

[0065]
Under the conditions
[0066]
[Expression 16]

[0067]
Use a weight that minimizes.
[0068]
Since this is the minimization of the sum of squares with the condition of the linear expression for wn, the solution can be easily obtained by solving the simultaneous linear equations by a known calculation method (Larange's undetermined multiplier method).
[0069]
<Weights for minimizing “folding components and noise”>
Considering the balance with the noise component, the aliasing component is also a kind of noise and does not necessarily need to be zero. A weight for reducing the aliasing component and the noise component is also effective.
[0070]
This does not make equation (26) a binding condition, but the square value of its absolute value as in the following equation:
[0071]
[Expression 17]

[0072]
Is regarded as a noise evaluation value equivalent to Equation (30), an evaluation function that is an increase function is taken for each noise evaluation value, and a weight that minimizes the evaluation function is used.
[0073]
Under the condition of Σwn = 1, for example, the weighted sum of Expression (32) and Expression (30) as the following expression is used as the evaluation function,
[0074]
[Formula 18]

[0075]
Use weights to minimize this. In this case, cancellation of the aliasing component has priority when α≈0, and noise reduction of the fundamental wave has priority when α≈1. Equation (33) is also the minimization of the sum of squares for wn, and a solution can be easily obtained by solving simultaneous linear equations. In this case, two or more sets of input data are required.
[0076]
<Sampling position estimation>
In the above description, the sampling position is considered to be known for each input sampling signal. However, when handling data captured while moving a digital camera, the sampling position cannot be known in advance, and it is necessary to estimate each sampling position from given digital data. The estimation method will be described below (in addition, the algorithm used here is “the distribution measurement system of velocity vectors using the spatio-temporal differential calculation method of the image”) (Proceedings of the Society of Instrument and Control Engineers) Vol.22, No.12, pp 88-94 (December 1986)).
[0077]
The deviation of the sampling position of the two digital data can be considered as the speed of the moving object. Therefore, a one-dimensional continuous function f (x, y, t) that translates at a speed (s, t) is considered as an example. If the change is small, as a first order approximation of the Taylor expansion,
[0078]
[Equation 19]

[0079]
Holds at each point. From the ∂f / ∂x, ∂f / ∂y and ∂f / ∂t in this equation, the change rate (s, t) can be obtained at each point.
[0080]
However, in the above formula,
(1) Deviation from first-order approximation
(2) Influence of noise
Therefore, the accuracy is poor. Therefore, a value obtained by integrating the square of the left side of the equation (34) within a certain area as an evaluation function is set as an evaluation function, and this evaluation function takes a minimum value (s, t). By calculating, an accurate value can be obtained.
[0081]
[Expression 20]

[0082]
Thus, if the relative position between two data is known, the relative position of all the data can be estimated by obtaining the relative position with respect to all the other data with one data as a reference.
[0083]
<Example>
FIG. 1 is a block diagram of an embodiment of the present invention. In this embodiment, based on the above-described principle of the present invention, two-dimensional digital data obtained by sampling the same object at a plurality of different sampling positions is input, and the high-frequency component of the original signal before being sampled is restored. A signal processing method for obtaining high-resolution digital data will be described. The function of each part will be described below.
[0084]
(Input digital data)
The input digital data 1 used in the present invention is two-dimensional digital data obtained by sampling the same object at the same sampling interval while shifting the position. For example, using a digital camera, digitized multiple data shot while moving the same object little by little, and multiple consecutive frames of a scene shot while moving the same object slowly with a video camera Data can be used. Further, by examining the MTF of the optical system of the input device such as a camera and the aperture characteristics of the sensor in advance, the band of the original signal before being sampled can be known. Here, an example for restoring the original signal band-limited to twice the Nyquist frequency is shown.
[0085]
Hereinafter, the input digital data are referred to as in0, in1,. The nth set of input digital data (in _n ) In i-th row and j-th column (pixel) _n Called [i, j].
[0086]
(Reference data buffer)
The reference data buffer 2 is a buffer for storing the 0th set of data (in0) to be a reference position of the sampling position. The subsequent position difference estimation unit inputs two sets of digital data and estimates the difference between the sampling positions. By using the 0th set of data (in0) as reference data and estimating the sampling position difference between in0 and in1 and in0 and in2, all the input data A relative sampling position can be estimated.
[0087]
(Position difference estimation unit)
The position difference estimation unit 3 estimates each input data (in1, in2,...) By estimating a sampling position difference from the reference data (in0) stored in the reference data buffer. Is a means for estimating the relative sampling position of.
[0088]
The sampling position of the first set of data is output as (p1, q1), the sampling position of the second set of data is output as (p2, q2), etc. (the sampling position of the reference 0th set of data is ( 0,0)). Such estimation can be realized by substituting the differential and integral of the equation (35) with the difference and addition of digital data.
[0089]
That is, in0 and in _n To estimate the position difference (pi, qi) of
[0090]
[Expression 21]

[0091]
Since the displacement of the sampling position is measured instead of the target speed (s, t), the sign is changed,
[0092]
[Expression 22]

[0093]
Is calculated as (pn, qn). This can be solved by the linear least square method.
[0094]
(Weight calculation part)
Based on the sampling positions ((0, 0), (p1, q1)...) Obtained by the position difference estimation unit 3, the weight calculation unit 4 obtains weights used in the subsequent product-sum unit. Means. .., W1,... Are output as weights for the 0th, 1.
[0095]
Specifically, here, since the original signal has a band twice the Nyquist frequency, nine sets of input digital data are used as described in the explanation of the principle, and (k, l = -1, 0, 1). As a result, an operation for solving the equations (22) and (26) is performed. As another method, 10 or more sets of input digital data can be used, and weights that minimize Expression (31) can be obtained under the constraints of Expressions (22) and (26). Further, as another method, two or more sets of input digital data can be used to obtain a weight that minimizes the expression (33) under the constraint condition of the expression (22). When the band of the original signal is different from this, a weight wn that cancels all unnecessary (k, l) harmonics is obtained.
[0096]
(Resolution magnification specification part)
The final output digital data in this embodiment is high-resolution data having a larger number of pixels than the input digital data. The resolution magnification specifying unit 5 is means for specifying the resolution of the output digital data as the resolution magnification of the input digital data. Here, it is assumed that the resolution is quadrupled. That is, the output digital data is four times as many as the input digital data and the sampling period is 1/4.
[0097]
The user may specify the resolution magnification. However, in the present embodiment, since up to twice the frequency component of the Nyquist frequency is restored, in order not to include aliasing distortion in the output data, the magnification must be twice or more.
[0098]
(Broadband LPF processing and high resolution processing unit)
This portion corresponds to the portion described in the principle of the present invention, “The signal obtained by adding all of the fundamental wave and a part of the adjacent harmonics is extracted by the LPF having a transmission band of [−1, 1]”.
[0099]
Input digital data is in _n [I, j] is a data string, and this is considered as a continuous function (delta pulse string) with respect to the position (x, y) as in the following expression. It corresponds to y).
[0100]
[Expression 23]

[0101]
This signal (dn (x, y)) is multiplied by the wideband LPF of equation (10).
[0102]
The ideal LPF process is in frequency space
[0103]
[Expression 24]

[0104]
In real space,
[0105]
[Expression 25]

[0106]
It is a convolution with. When a signal obtained by LPF processing of dn (x, y) is en (x, y),
[0107]
[Equation 26]

[0108]
It is.
[0109]
Furthermore, the input digital data (in _n [I, j])
[0110]
[Expression 27]

[0111]
Can be written.
[0112]
A signal (en (x, y)) obtained by this LPF processing is a continuous signal, and a signal obtained by sampling this is used as an output digital signal of the block 6. Therefore, it is only necessary to calculate the equation (38) for the desired (sampled) position (x, y) in the continuous signal.
[0113]
The digital signals obtained here are added by the product-sum unit 8 at the subsequent stage. Addition is to add pixel values of the same position to the original signal. Therefore, the position where en is sampled is preferably a position corresponding to the same position on the original image on each input image having a different sampling position.
[0114]
Here, the pixel position of the high resolution data mid0 [i, j] of the reference data in0 [i, j] is used as a reference, and the other input data corresponds to the pixel position of the high resolution standard data. The pixel value at the position on the image is obtained.
[0115]
FIG. 8 shows the positional relationship of high resolution pixels. In the figure, large white circles indicate pixels of reference input data sampled at the sampling interval 1. A large black circle ● indicates a pixel position of the nth set of input data sampled at the sampling interval 1 at a position shifted by (pn, qn) with respect to the reference input data.
[0116]
First, the high resolution processing unit simply sets the sampling interval to 1/4 with respect to the reference input data, and obtains a pixel value at a position (i / 4, j / 4) (i and j are integers). . In FIG. 8, small white circles indicate pixel positions obtained by increasing the resolution of the reference input data.
[0117]
Nth set of input data (in _n The resolution-enhanced data is created as follows.
That is, the pixel value corresponding to the pixel position (small white circle in FIG. 8) of the high-resolution reference data on the original signal using the pixel value (black circle) of the n-th set of input data (for example, interpolating). Are obtained as n-th set of high resolution data. FIG. 9 shows pixel positions of the n-th set of high resolution data created by, for example, interpolating the pixel values of the n-th set of input data. The position of the small black circle in the figure is the same position as the pixel position (small white circle) obtained by increasing the resolution of the reference input data. Since the sampling position of the n-th set of input data is shifted by (pn, qn), ((i / 4) −pn, (j / 4) −qn) is a corresponding position.
[0118]
In the previously proposed method, the pixel position of the resolution-enhanced data is determined relative to each input data, and the position is shifted and added when obtaining the weighted sum thereof. In this method, the unit of shifting the pixel position for addition is the pixel interval of the high resolution data. For example, when the magnification of the resolution-enhanced data is set to 4 times, it can be shifted only by 1/4 pixel unit of the input data. Therefore, when the magnification is small, the alignment accuracy is poor, and as a result, the accuracy of the restored signal may be deteriorated. On the other hand, the method of the present invention is based on the n-th set of input data (that is, pixel values) and the sampling position (pn, qn), and the same position as the pixel position of the reference high-resolution data. Therefore, the pixel position of the high-resolution data can be determined with high accuracy regardless of the magnification for high resolution.
[0119]
The digital data that is output from the “wideband LPF processing and high resolution processing unit” obtained in this way is called midn [i, j]. Here, the new sampling interval is 1/4, so midn [i, j] is
[0120]
[Expression 28]

[0121]
If a table called sinc2 (i, j) (actually an approximation of it within a finite range) is prepared in advance, the input digital signal in _n Can be calculated by multiplying the products.
[0122]
The LPF-processed signal en (x, y) is a signal including a frequency component higher than the Nyquist frequency of the original signal. In order to express it effectively, that is, to restore the high frequency component of the original signal, the sampling interval for the LPF-processed signal en (x, y) must be narrower than the sampling interval of the input signal. I must. If the sampling interval for the signal en (x, y) is not narrower than the sampling interval of the input signal, aliasing distortion occurs again, which is meaningless. Therefore, it is set to 1/4 here. Further, when the band of the original signal is different from this, an LPF that transmits the band necessary for restoring the original signal is used.
[0123]
(Intermediate data buffer)
The intermediate data buffer 7 is a buffer for storing intermediate digital data (midn [i, j]) obtained by the “broadband LPF processing and high resolution processing unit” in the previous stage.
[0124]
(Product sum)
The weighted sum calculation in the product-sum unit 7 corresponds to the part “takes the weighted sum of Expression (12) and cancels the harmonic component” described in the principle of the present invention. In the configuration of this embodiment, processing is performed in real space. In the real space, each intermediate data may be added after weighting wn.
[0125]
Since the intermediate data buffer 7 stores data corresponding to the same position on the original signal, the pixel values (mid0, mid1...) Of the respective buffers are multiplied by weights w0, w1. Just be happy.
[0126]
The final output data (out [i, j]) is
[0127]
[Expression 29]

[0128]
It becomes.
[0129]
【The invention's effect】
As explained above, The present invention According to the present invention, for the n sets of two-dimensional digital data obtained by sampling the same signal n times at the same sampling interval while changing the sampling position, a filter that transmits all necessary bands of the original signal, By taking a weighted sum using weights that cancel all unnecessary harmonics that have passed through the filter, it is possible to cancel all unnecessary aliasing distortions due to sampling and restore all high-frequency components of the necessary original signal. This makes it possible to obtain high-resolution data with less blur and distortion.
[0130]
In addition, the present invention According to this, since the aliasing component is completely canceled with respect to data that does not include noise, the high-frequency component of the original signal can be restored with high accuracy.
[0131]
In addition, the present invention According to the above, the aliasing component and the noise can be simultaneously reduced even for the data including the noise.
[0132]
In addition, the present invention According to this, it is possible to accurately cancel the aliasing component and restore the high-frequency component of the original signal, regardless of the magnification for increasing the resolution.
[Brief description of the drawings]
FIG. 1 is a block diagram of an embodiment of the present invention.
FIGS. 2A, 2B, and 2C are diagrams illustrating a sampling process in a real space. FIG.
FIGS. 3A, 3B, and 3C are diagrams illustrating a sampling process in a frequency space.
FIG. 4 shows an original signal band limited to twice the Nyquist frequency.
FIG. 5 shows a signal in which aliasing has occurred due to sampling.
FIG. 6 shows an LPF processed signal.
FIG. 7 shows a signal in which triple folding occurs due to sampling.
FIG. 8 shows a positional relationship of high resolution pixels.
FIG. 9 shows the positions of n-th set of high resolution pixels.
[Explanation of symbols]
1 input digital data
2 Reference data buffer
3 Position difference estimation unit
4 Weight calculation part
5 Output resolution specification part
6 Wideband LPF processing and high resolution processing section
7 Intermediate data buffer
8 Sum of products
9 Output digital data

Claims

同一の信号について、標本化位置を変えて同一の標本化間隔でｎ回の標本化によって得られるｎ組の２次元デジタルデータに対して、標本化周波数の１／２以上の周波数帯域を含む、必要な原信号の帯域を全て透過するローパスフィルタ処理を行い、該ローパスフィルタ処理によって得られたｎ組の連続信号をそれぞれ前記標本化間隔より狭い間隔で標本化した画素値を高解像度化データとして求める高解像度化処理手段と、前記各デジタルデータの標本化位置に応じて、前記透過帯域に含まれる不要な折り返し成分を打ち消すような重みを算出する手段と、対応する位置における各高解像度化データに該重みをつけて加重和をとり出力する手段とを備えたことを特徴とする２次元デジタル信号処理装置。For n sets of two-dimensional digital data obtained by sampling n times at the same sampling interval at the same sampling interval with respect to the same signal, including a frequency band of 1/2 or more of the sampling frequency. Low-pass filter processing that transmits all necessary bands of the original signal is performed, and pixel values obtained by sampling n sets of continuous signals obtained by the low-pass filter processing at intervals smaller than the sampling interval are used as high-resolution data. high and resolution processing means, in accordance with the sampling position of each digital data, means for calculating a weight for canceling the unwanted aliasing components included in the transmission band, each high-resolution data in the corresponding position to seek A two-dimensional digital signal processing apparatus comprising: a means for adding the weight to the sum and outputting the weighted sum.