DIYウィンドウ機能

デジタル信号処理では、ウィンドウ関数は信号を時間内に制限するために広く使用されており、その名前は、ハン、ハミング、ブラックマン、ハリスなど、何らかの方法で離散フーリエ変換に出くわしたすべての人によく知られています。しかし、それらは十分であり、何か新しいものを思い付くことが可能であり、それに何か意味がありますか?



この記事では、WolframMathematicaを使用して新しいプロパティを持つウィンドウ関数の出力を見ていきます。また、読者は、議論中の問題のコンテキストでのデジタル信号処理の一般的な理解があり、少なくともWikipedia記事に精通していることを前提としています







前書き



歴史的に、最初のウィンドウ関数は、サイドローブの振幅を減らすことによってスペクトル特性を改善しようとする過程で登場しました-信号にウィンドウ関数を掛けると、スペクトルが複雑になり、分析にいくつかの制限が生じます。 20世紀の50年代の当時、計算能力はシンボリック計算を簡単かつ自然に操作することを可能にせず、最適なパラメーターの選択は非常に困難でした。これが、関数が異なる名前で呼ばれる理由の1つです。つまり、それらのスペクトル特性は、かなり長い間、さまざまな研究者によって研究されてきました。これの副作用は、名前付きウィンドウ関数の既存のセットが一種の「標準セット」として認識され、それ以外では何も見つからないことです。同時に、これらのウィンドウの名前にはプロパティに関する情報が含まれておらず、個別の調査と記憶を示唆しています。



分類



すべてのウィンドウ関数は、一部の関数に長方形の関数を乗算することによって取得されます。これにより、時間の制約が保証されます。したがって、まず第一に、それらはそれらが使用する機能に従って分類することができます:



  1. コサインの合計。それらのスペクトルは簡単に計算可能であり、重み付けされたsinc関数の合計であるという事実により、最も広範なクラスです。これには、ハン、ハミング、ブラックマン、ハリスが含まれます...
  2. ピースワイズ連続多項式。それらは、単純な関数(たとえば、長方形や三角形)の畳み込みの結果として取得されます。同時に、それらのスペクトルは乗算され、それらの発見も特に困難ではありません。これらには、Dirichlet、Triangular、Parzen、Welchなどが含まれます...
  3. 他のすべて-指数、ガウス、sincなどを使用して、特定のスペクトル特性よりも本質的にイデオロギー的である特定の関数の選択。


それらは、エッジのプロパティで分割することもできます。



  1. エッジにギャップはありません-値はゼロに等しくなります。休憩はディリクレ、ハミング、ブラックマン-ハリスです。持っていない-三角、ハン、ヌタール
  2. 一次および他の派生物の不連続性の欠如


さらに2つのウィンドウ機能が別々に目立ちます。



  1. カイザー。メインの花びらの幅を設定できます。
  2. Dolph-Chebyshev、そのすべてのサイドローブは与えられた振幅に等しい。


それらの両方は、値と導関数の両方のエッジに不連続性があり、それらの計算にはいくつかの困難が伴います-カイザー関数は特別な関数(ベッセル)によって計算され、ドルフ-チェビシェフ関数は周波数領域で決定され、逆離散フーリエ変換によって計算されます それらの分析的抗誘導体を見つけることも特に困難です。



次の簡単なコードを使用して、ウィンドウ関数の中断を調べることができます。
wins = {DirichletWindow, HammingWindow, BlackmanWindow, 
   BartlettHannWindow, BartlettWindow, BlackmanHarrisWindow, 
   BlackmanNuttallWindow, BohmanWindow, ExactBlackmanWindow, 
   FlatTopWindow, KaiserBesselWindow, LanczosWindow, NuttallWindow};
Table[{f, 
  Table[Limit[D[f[x], {x, k}], x -> 1/2, 
    Direction -> "FromBelow"], {k, 0, 4}]}, {f, wins}]</code>
<img src="https://habrastorage.org/webt/rm/wq/8f/rmwq8fjkjxjcox-czh8zevkbd0m.png" />

<code>Plot[Evaluate[Through[wins[x]]], {x, -1, 1}, PlotRange -> {-0.25, 1}, 
 GridLines -> {{-0.5, -0.25, 0.25, 0.5, 1}}, AspectRatio -> 5/8, 
 PlotLegends -> "Expressions"]








リバースエンジニアリング



Blackman関数とNuttel関数の定義を見てみましょう。



BlackmanWindow[x] // FunctionExpand




{150(25cos(2πx)+4cos(4πx)+21)12x120|x|>12



NuttallWindow[x] // FunctionExpand




{121849cos(2πx)+36058cos(4πx)+3151cos(6πx)+8894225000012x120|x|>12



それらは、偶数の周波数(ゼロから開始)といくつかの係数を持つ3つと4つの余弦波の合計です。これらの係数はどこから来たのですか?それらが手動で選択された可能性は低く、少なくともそれぞれ個別に-結局のところ、スペクトル特性に関係なくウィンドウが準拠しなければならない境界条件があります-少なくとも中央の1つに等しいです。



ブラックマン機能を自分たちで「発明」してみましょう。これを行うために、まだ未知の係数を持つ関数を定義します



f = Function[x, a + b Cos[2 Pi x] + c Cos[4 Pi x]];


そして、境界条件を定義する方程式のシステムを構成します-中心で1に等しく、端でゼロに、そして端で導関数のゼロに等しく。



ss = Solve[{
   f[0] == 1,
   f[1/2] == 0,
   f'[1/2] == 0
   }, {a, b, c}]








{{b12,c12a}}



解決策は見つかりましたが、1つではなく、多くあります。係数aに応じて、Wolframは丁寧に警告しました。ここで、見つかった解から、未知の係数に応じて新しい関数を定義します。a = 0.42の場合、ブラックマン関数が得られる ことは簡単にわかります。しかし、なぜ正確に0.42なのですか? この質問に答えるには、そのスペクトルをプロットする必要があります。分析的フーリエ変換は最も簡単な作業ではありませんが、Wolframもそれを処理できるため、雑用から解放されます。



hx = Function[{x, a},

Piecewise[{{f[x] /. ss[[1]], -1/2 < x < 1/2}}, 0] // Evaluate]
















hw = Function[{w, a}, 
  FourierTransform[hx[x, a], x, w] // #/Limit[#, w -> 0] & // 
    Simplify // Evaluate]






注意
, FourierTransform , — , , — . w. — w, , .







この形式では、対数スケールでスペクトルを分析する方が便利なため、関数のグラフはまだあまり有益ではありません。デシベルに適切な変換を追加することによって20log10(|x|)、ウィンドウ関数のスペクトルの同じ通常のグラフを取得します。



コード
Manipulate[
 Plot[{
      hw[w, a],
      hw[w, 0.42]
      } // Abs // 20 Log[10, #] & // Evaluate
  , {w, -111, 111}, PlotRange -> {-120, 0}, GridLines -> Automatic, 
  PlotStyle -> {Default, Thin}, PlotPoints -> 50]
 , {{a, 0.409}, 0.3, 0.5}]








パラメータを変更する過程で、同様のことが観察さ







れます。パラメータaに応じて、サイドローブのレベルが変化し、a = 0.42では、ほぼ最小で均一です。では、我々はサイドローブの可能な最小レベルを意味し、「より良い」ことであれば= 0.409、我々は若干良い結果を得ることができます。



注意
, , .



開発



明らかに、そのようなマイクロ最適化は努力する価値がまったくありませんでした。初期データを変更せずに、このウィンドウのプロパティを定性的に改善することは可能ですか?



先ほど、合計が1になるウィンドウ関数の出力を確認しました。これにより、元の信号を復元できます。この手法をウィンドウで試して、スペクトルにどのように影響するかを見てみましょう。



-1/2から1/2の範囲での引数の変化の境界を考慮して、オーバーラップするウィンドウを構築するための補助関数を定義しましょう。







抗誘導体を見つけ、座標の中心にシフトして1にスケーリングし ます











グラフに表示します。







ご覧のとおり、ここでもWolframが独自に実行し、抗誘導体の断片的な連続定義を手動で指定する必要はありませんでした。これで、ウィンドウのビューは変数aだけでなく、オーバーラップの程度にも依存します。ウィンドウが大きくなると、派生物の形式になる傾向があります。



コード
Manipulate[
 Plot[
  OverlapShape[hsx, x, a, t]
  , {x, -1, 1}, PlotRange -> All, GridLines -> Automatic]
 , {{a, 0.404}, 0.4, 0.5}, {{t, 4}, 3, 11}]






そして最後の仕上げは、パラメータaの最適値を決定するために、スペクトルの分析関数を見つけることです。



ここで、前回のように変換を直接計算しようとすると、Wolframを深く考えさせます。このプロセスを高速化するには、いくつかの方法があります



。-分析変換の代わりに離散変換を使用する-最も単純な変換。しかし、実際の数学者は、純粋に分析的な解決策を見つけることができる数値的方法を適用しません-そして私たちは(まだ)しません。さらに、明らかな制限があります(解像度とスペクトルの制限に関して)。



-オーバーラップパラメータとして特定の値を渡し、結果を見て、パターンを確認して一般化してみてください。かなり実用的なオプションですが、追加の精神的および創造的な努力が必要です。最後の手段として残しましょう。



-周波数領域で直接計算します。これは、フーリエ変換の次の特性により可能になります



。1)線形、つまり画像の合計は、合計の画像と同じです。

2)時間領域での積分は、周波数領域でのiw除算することと同じです。

3)時間変位は、複雑な正弦波による変調と同等です。

ウィキペディアまたはここで詳細)。



したがって、任意のウィンドウ関数のスペクトルhwを取得すると、tがオーバーラップする関数のスペクトルhwoを取得 できます











これ







で、ダイナミクスでスペクトルがどのように変化するかを確認できます。ここで、オーバーラップパラメータを変更すると、ウィンドウスペクトルにわずかに異なる方法で影響します。 :







パラメータを少し試してみると、「多ければ多いほど良くない」ことに気づきやすく、その関数の最適なオーバーラップの程度は4の領域にあります。具体的には、t = 4およびa= 0.404 -80dBを超えないサイドローブレベルが得られます。これは非常に良い結果です。特に、伝統的に合計ウィンドウに使用される上げられたコサイン関数が約-30dBのローブレベルを与えることを考えると。明示的に書き出すと、新しいウィンドウ関数は次のようになります。



コード
OverlapShape[hsx, x, 404/1000, 4] // PiecewiseExpand // FullSimplify




{404π(12x)+36sin(13(16πx+π))+375sin(13(π8πx))606π14<x12404(2πx+π)+36sin(23π(8x+1))+375sin(13(8πx+π))606π12<x143(125cos(8πx3)+12cos(16πx3))202π+1314<x140|x|>12



さらなる開発



サイドローブレベルをさらに下げるために他に何ができますか?偶数ではなく奇数の周波数でコサインを取ることができます(ここでは、コンパクトにするために、方程式系の解は関数の定義に直接埋め込まれています):



コード
hx1 = Function[{x, a},
  a Cos[1 Pi #] + b Cos[3 Pi #] + c Cos[5 Pi #] & // #[x] /. Solve[{
          #[0] == 1,
          #[1/2] == 0,
          #'[1/2] == 0,
          #''[1/2] == 0
          }, {a, b, c}][[1]] & // # UnitBox[x] & // Evaluate]




(acos(πx)+(58a2)cos(3πx)+(38a2)cos(5πx))Π(x)



パラメータa = 0.6628およびオーバーラップレベル4.5と統合した後 -90dBの減衰を取得します。







明示的に:

{327sin(17π(45x+2))+24855sin(17π(19x))+3670cos(114π(54x+1))+2151243024518<x1224855sin(17π(9x+1))+3670sin(37π(9x+1))+327sin(57π(9x+1))+215124302412<x5183670sin(37π(9x1))+24855sin(17π(9x+1))+327sin(57π(9x+1))43024++3670sin(37π(9x+1))+327sin(17π(45x+2))+24855sin(17π(19x))43024518<x5180|x|>12



コサインをもう1つ追加して、ゼロ導関数の数を増やすことができます。



コード
hx2 = Function[{x, a},
  a Cos[1 Pi #] + b Cos[3 Pi #] + c Cos[5 Pi #] + 
       d Cos[7 Pi #] & // #[x] /. Solve[{
          #[0] == 1,
          #[1/2] == 0,
          #'[1/2] == 0,
          #''[1/2] == 0,
          #'''[1/2] == 0,
          #''''[1/2] == 0
          }, {b, c, d}][[1]] & // # UnitBox[x] & // Evaluate]




(acos(πx)+180(3516a)cos(3πx)+180(3548a)cos(5πx)+140(58a)cos(7πx))Π(x)



パラメータa = 0.5862およびオーバーラップレベル6.4と統合した後、-110 dBの抑制を取得します。







明示的に:

{3077550sin(154π(64x5))90069sin(554π(64x5))5820sin(754π(64x5))560455cos(19π(732x))+260134452026881132<x12560455sin(118π(64x+5))+3077550sin(154π(64x+5))+90069sin(554π(64x+5))+5820sin(754π(64x+5))+2601344520268812<x11323077550sin(154π(64x5))90069sin(554π(64x5))5820sin(754π(64x5))560455cos(19π(732x))5202688++560455sin(118π(64x+5))+3077550sin(154π(64x+5))+90069sin(554π(64x+5))+5820sin(754π(64x+5))52026881132<x11320|x|>12



サイドローブのレベルをさらに大幅に下げるには、フーリエ変換を2乗して、サイドローブのレベルを一度に2分の1に減らします。これにより、手動で選択するためのパラメーターの数の増加を取り除くことができますが、時間領域での畳み込みの計算が複雑になります。



コード
hx3 = Function[{x, a}, 
  Convolve[hx1[x, a], hx1[x, a], x, y] /. y -> 2 x
 // FullSimplify // Evaluate]








そして、ここではすでに160 dBを超える抑制を達成することが可能







です。その印象的なサイズのため、式を明示的に示しません。



ハイパージオメトリックウィンドウ関数



関数の境界で必要な数のゼロを確保するために、方程式系の解を検索してから積分しました。これはあまり便利ではありません。毎回方程式の数を変更する必要があります。たぶん、もっとシンプルで美しい解決策がありますか?有る!そして、ハイパージオメトリック関数はこれに役立ちます。



ハイパージオメトリック関数について簡単に
, . — , . — .



私たちのソリューションは次のようになります。

2zΓ(n+32)2F1(12,n2;32;z2)πΓ(n2+1)+a(z+bz3)(1z2)n/2

どこ z=sin(πx2)



それは2つの部分の合計で構成されています:関数のアンチデリバティブ (1z2)n2振幅が点(1,1)に正規化され、自由なパラメータを持つ関数に同じ値が乗算されます (1z2)n2必要な数のゼロ導関数を維持するために:







調整精度とサイドローブの分布への影響は、形状を修正するために選択した機能によって異なります。個別の調査が必要なオプションがあります。この場合、関数a(z+bz3)完全に受け入れ可能なオプションを提供します-2つのパラメーターにより、より正確な設定を行うことができます。



数式自体は、偶数nの場合は多項式に簡略化され、奇数nの場合はアークサインと平方根を持つ多項式の合計に簡略化されるという点で興味深い(そして便利です)。これにより、最終的な数式がよりコンパクトになり、計算が簡単になります。



離散フーリエ変換を介してここでパラメータを選択することは、すでに簡単(かつ高速)です。また、3つのパラメーターを操作するには、オーバーラップ関数の追加の定義が必要です。



















例として、図のパラメーターを置き換えた後、関数は次のように簡略化されます。



コード


288z11+2315z97380z7+12330z511940z3+8163z3200





逆カイザーウィンドウ



ウィンドウを1つにまとめる作業に値しない場合、最適なウィンドウはカイザーウィンドウです。これにより、メインローブの幅をスムーズに調整できます。ただし、これには欠点があります。ベッセル関数で表されるため、数学パッケージの外部で読み取るのはやや困難です。もちろん、ベッセル関数を個別に実装することも、基本関数を使用して近似値を探すこともできます。そして突然、カイザーウィンドウのスペクトルの関数(時間制限)を使用すると、さらにわずかに良い結果を得ることができることがわかりました-



sinh(k14x2)sinh(k)14x2Π(x)

注意
±12 , ksinh(k).


そのような狡猾な決定に対する支払いは、そのスペクトルを純粋に分析的に表現することは不可能であることが判明しました-しかし、実際には、いずれにせよ、離散フーリエ変換を使用して離散形式でデータを操作および分析するため、これはそれほど重要ではありません。下のグラフでは、スペクトルを比較できます。ここで、赤は元のカイザーウィンドウ、緑はその近似値です。



コード








メインの花びらの幅が同じであるため、サイドローブのレベルはわずかに低くなります。また、使いやすさのために、パラメータkの概算値を使用してテーブルを作成し、必要なレベルのサイドローブ抑制を取得できます

-60 k=8.8-90 k=11.36-120 k=15.18-150 k=18.88



逆カイザーウィンドウのアンチデリバティブの近似値を見つけることは、別の興味深い問題のままです。これまでのところ、それを表現することは、かなりゆっくりと収束するパワーシリーズによってのみ可能でした。



sinh(k14x2)sinh(k)14x2dx=n=1212nkn1x2n1(k(1)nKn12(k)+kKn12(k))π(2n1)sinh(k)Γ(n)



おそらく数学の専門家はコメントでより良い解決策を提案することができます。



結論



デジタル信号処理は完全に形成された分野であるという事実にもかかわらず、開発の余地があり、新しいものを見つける余地があります。少なくとも既存の特別な科学文献では、1つのユニットにまとめられたウィンドウ機能の構築は考慮されていません。また、最新のコンピューター代数システムの機能により、蓄積された知識を質的に新しいレベルで再考することができます。



すべての計算を含む元のWolframMathematicaドキュメントは、GitHubで入手できます



All Articles