この投稿は、 スパース時間周波数表現を介したアダプティブデータ分析の記事で説明されている比較的新しい信号処理方法と、個人的に私を混乱させた小さな間違いに関するものです。 この記事は、2011年にカリフォルニア工科大学の応用数学の教授であるThomas I. HoweとShi Zotzzyanによって出版されました。おそらく、あなたがこれを読んでいる頃にはすでに修正されています。
この記事に出くわしたのは、非線形信号と非定常信号の時間周波数解析のさまざまな方法を探していたときです。私の場合は、人間の血管内の均一な血液要素を動かすことによる超音波信号です。 この分析の本質は、信号の特性の変化を追跡することです。つまり、信号の周波数の時間依存性を知りたいのです。 広範囲に及ぶ方法を除いて-スペクトルおよびウェーブレット分析、EMD（経験的モードへの分解）および正弦波モデリングなどの方法が見つかりました。これについては後で説明します。
経験的モードの方法は非常に簡単に使用できますが、結果の妥当性に関して特に開発されたわけではありません。 Thomas HouとShi Zuoqiangは数学的装置の開発をさらに進め、正弦波信号モデリングの独自の方法を提案しました。 彼のアイデアは、滑らかな振幅を持つ高調波への信号のスパース分解です。 期待される結果は上の写真にあります。 この場合、関数f（t）= 6t + cos（8πt）+ 0.5 cos（40πt）によって得られた信号が分解されました。 もちろん、信号の分解は一意ではないため、成分高調波の最小基準が導入され、次のように問題が形成されました。
問題P0の解決策は、振幅aと位相関数thetaの高調波であり、その微分係数は信号周波数に対応します。 著者は、この問題を再帰的に解決することを提案しています。1つの高調波成分を抽出し、残り（中央値）を取得します-それを処理し、次の高調波を抽出し、残りを取得します。 中央値は信号よりもはるかに滑らかである必要があるため、高調波の検出は次の問題を解決するために削減されます。
ここで、a1 * cos（theta1（t））は得られた高調波、a0は中央値です。 振幅の滑らかさは、合計変動によって決まります。
関数gの滑らかさが大きいほど、その（n + 1）導関数がゼロに近くなり、それに応じて、n次の全体的な変動が小さくなります。 これを言うこともできます：補間できる多項式の次数が小さいほど、関数は滑らかになります。 最小のゼロ次変動を探すのが最善であると仮定するのが論理的ですが、この場合、区分的に一定の関数（階段効果）を取得し、必要な高次の情報（曲率など）を失うリスクがあります。 この記事の著者は、3次変動を減らし、4次導関数をゼロにすることを提案し、それによって3次スプラインによる近似に似たものを作成します。 その結果、最小TV ^ 3（a0）を探しています。 しかし、a1が十分に滑らかであることも重要であるため、最良の場合、最小TV ^ 3（a0）+ minλ* TV ^ 3（a1）を見つける必要があります。ここで、λはラグランジュ乗数で、この場合は重み係数の役割。 簡単にするために、著者はλ= 1のままにして、次の問題に到達しました。
アルゴリズム
この問題を解決するために、反復アルゴリズムが提案されています。 この条件では、1つ以上の高調波が追加されます-b1 * sin（theta1（t））、その滑らかさも目標になります。 条件関数cos（theta ^ 0）= f（t）を満たす増加関数になるように位相関数theta ^ 0を初期化し、メインの反復を開始します。
最初のステップでは、振幅の一般的な変動を最小限に抑えます。 次に-上記の式に従ってthetaを更新し、どれだけ変化したかを確認します。 アルゴリズムの終わりまでに、振幅b1はそれぞれ無視できるようになり、artcg（b1 / a1）は小さくなり、シータはある値に収束すると想定されます。 これがフェーズ関数になります。 シータ^ nを受け取った後、最初のハーモニックa1 * cos（theta1（t））を区別します。ここで、シータ1 =シータ^ nです。 次に、このアルゴリズムを剰余に対して実行します。つまり、値f（t）= a0 + b1 * sin（theta1（t））を使用します。 次の高調波a2 * cos（theta2（t））を選択し、残りを続行します。 そして、残りの標準が無視できるようになるまで続きます。 私は最初の規範を取りました。 その定義を覚えていない人のために：
総変動の最小値を見つける問題を適用するために、著者はL1最小化の問題に還元し、微分方程式を差分方程式に置き換えました。 信号、または拡張したいNポイント（ノード）の順序セットがあります。つまり、ベクトルa0、a1 * cos（theta1（t））およびb1 * sin（theta1（t））の合計として存在します。 そして、再定式化された問題では、不幸なタイプミスが認められました。
a0 ^ n、a1 ^ n、b1 ^ nがグリッド関数であると仮定すると、それらに行列D ^ 4を掛けて4次導関数を近似します。 次に、取得した値のモジュールの合計を最小化して、対応する変動を減らします。 すべてが論理的で信able性がありますが、このアルゴリズムを実装した結果、私は結果を見るのを恐れました-それは完全に混oticとしたものであることが判明しました。 長い間コード内のエラーを見つけようとして、私はついにそのエラーが記事にあることに気付きました。
実際には、Phiにxを掛けることにより、対応する点での要素が導関数a0 ^ n、a1 ^ n、b1 ^ nの合計であるベクトルが得られます。 つまり、上記の3つのバリエーションの合計ではなく、（a0 ^ n + a1 ^ n + b1 ^ n）の合計バリエーションを最小化します。
エラーはエラーであり、発見され、修正され、数か月間忘れられました。 そして、5月の休日に、私はまだ記事min。| Phi * x ||で置き換える提案をして、教授に手紙を書くために引っ張りました。 最小|| D ^ 4 * a0 ^ n || + || D ^ 4 * a1 ^ n || + || D ^ 4 * b1 ^ n ||。 驚いたことに、1時間後、ハウ教授から「メールをありがとう。 私たちはあなたのメモを勉強し、すぐにあなたに連絡します。そして、簡単な「トム」の署名で。
さらに6時間後、志教授から手紙が届いた。
親愛なるアレクサンドル、
メールありがとうございます。 ここでタイプミスをしました。
\ Phiはdiag（D4、D4、D4）でなければなりません。
アルゴリズム。 ありがとう。
最高
強強市
当然、私の改正と彼らが書きたかったことは同一です。 私が注意しなければならない唯一のことは、それがスパースであるという仮定の下でPhi行列を格納しない場合、つまりゼロ要素の束を格納する場合、メモリはすぐに終了する可能性があるということです。
L1最小化問題を解決する方法
アルゴリズムについては、教授は分割法Bregman反復法とL1正規化最小二乗法を提案しています。 ただし、これらの方法はかなり最近のものであり、おそらくすべてのパッケージが配線されているわけではありません。 怠け者（私のように）にとっては、最初のノルムの最小値を見つけるアルゴリズムを実装して、この問題がより一般的な線形計画問題にどのように還元されるかを知るのに役立ちます。 元のタスクは次のとおりです。
の対象
ここで、Aは行列NxM、Phiは行列KxM、xは次元Mのベクトルです。
これは、次の線形プログラミングタスクと同等です。
条件の下で
ここで、phi（i、j）は行列Phiの要素です。
正式には、移行は次のように実行できます。 K変数がベクトルxに追加され、それにより次元K + Mになります。
次元NxKの行列が行列Aに追加されます。
ベクトルc（LP問題の係数）は、M個のゼロとK個のユニットで構成されます。
行列Fは、行列Phiと単位行列Iで構成されます。次元F-（2 * K）x（M + K）：
したがって、正式なタスクは次のようになります。
条件の下で
結果
前述の例f（t）= 6t + cos（8πt）+ 0.5 cos（40πt）のアルゴリズムの結果を次に示します。 グラフでは、このアルゴリズムによって得られた値（青）と分析値（赤）および経験的モードに分解して得られた値（黒）と比較されます。
次の例は、この方法が信号の周波数変化を追跡するのにどのように役立つかを示しています。
他の例は、実際のデータで動作するアルゴリズムの結果と同様に記事に示されています。
まとめ
このアルゴリズムは非常に興味深いものであり、私の意見では、さらなる分析が必要です。 その明らかな欠点は、計算の複雑さ、ノイズに対する感度、および結果としてノイズの多いデータでの不安定な動作です。 著者らは、ローパスフィルタリングとフーリエ変換を使用して感度の問題を解決することを提案しています。 しかし、その別の時間についての詳細。
最後に、私はもう一度自慢し、Kaltekhaの教授が取るに足りないものではあるが間違いを起こし、まだ学士を持っていないロシアの学生がこれを指摘できることは明らかであると言いたいと思います。