調和振動

フーリエ変換に関するHabréの記事や、デジタル信号処理(DSP)などのあらゆる種類の美しいものに関する記事がいくつかありましたが、経験の浅いユーザーは、なぜこれが必要なのか、どこで、そして最も重要なのはそれを適用する方法です。





ノイズの周波数応答。



個人的には、これらの記事(たとえば、 この記事)を読んだ後、それが何であり、なぜ実生活で必要なのかは明らかになりませんでしたが、面白くて美しいものでした。

美しい写真を見るだけでなく、そのように言って、何がどのように機能するかを心から感じたいと思います。 そして、オーディオファイルの生成と処理の具体例を示します。 音を聞いてそのスペクトルを見て、なぜそうなのかを理解することが可能になります。

この記事は、複雑な変数の機能理論、DSP、その他の恐ろしいトピックを所有している人にとっては興味深いものではありません。 好奇心の強い、学童、学生、そして彼らに同情する人にとっては、より可能性が高い:)。





私はすぐに予約をします、私は数学者ではありません、そして多くのことを間違って言うことさえできます(個人的なメッセージによって訂正します)、そして私は現在のプロセスの私自身の経験と私自身の理解に基づいてこの記事を書いています。 準備ができたら、行きましょう。



材料についてのいくつかの言葉





数学の学校のコースを思い出すと、サイングラフをプロットするために円を使用しました。 一般に、回転運動は正弦波に変換できることがわかります(あらゆる調和振動のように)。 このプロセスの最良の実例はウィキペディアにあります。





調和振動



つまり 実際、サイングラフはベクトルの回転から取得されます。これは、式:



f(x) = A sin (ωt + φ),







ここで、Aはベクトルの長さ(振動の振幅)、φはゼロ時間でのベクトルの初期角度(位相)、ωは回転の角速度です。



ω= 2πf、ここでfはヘルツ単位の周波数です。



ご覧のように、信号の周波数、振幅、角度がわかっていれば、高調波信号を作成できます。



魔法は、絶対にあらゆる信号の表現がさまざまな正弦波の合計(多くの場合無限)として表現できることが判明したときに始まります。 言い換えれば、フーリエ級数の形で。

英語版ウィキペディアから例を挙げます。 たとえば、のこぎり波信号を取得します。





のこぎり



その量は、次の式で表されます。







続けて要約すると、最初にn = 1、次にn = 2などとなると、高調波正弦波信号が徐々にノコギリに変わる様子がわかります。







おそらく最も美しいのは、ネットワークのオープンスペースで見つけた1つのプログラムです。 正弦グラフは回転ベクトルの投影であると既に述べましたが、より複雑な信号はどうでしょうか? これは、奇妙なことに、回転ベクトルのセットの投影、またはそれらの合計であり、すべて次のようになります:





ベクトルは、のこぎりを描きます。



一般に、 リンクたどってパラメータを自分で試してみて、信号がどのように変化するかを確認することをお勧めします 。 私見理解のためのより視覚的なおもちゃ、私はまだ会っていない。



フーリエ変換と呼ばれる、特定の信号から周波数、振幅、初期位相(角度)を取得できる逆手順があることにも注意してください。





いくつかの既知の周期関数のフーリエ展開( ここから



これについては詳しく説明しませんが、これを生活にどのように適用できるかを示します。 参考文献のリストで、材料について詳しく読むことができる場所をお勧めします。



実践的な演習に進みます!





すべての学生が講義に座っている間、たとえばマタンで質問をしているように思えます。なぜこのナンセンスが必要なのですか? そして、原則として、予見可能な将来に答えを見つけられないことは、残念ながら、主題への関心を失います。 したがって、私はすぐにこの知識の実用的な応用を示します、そして、あなたはすでにあなた自身でこの知識を習得します:)。



その場でさらにすべてを実装します。 もちろん、Linuxの下ですべてを行いましたが、仕様は一切使用しませんでした;理論的には、プログラムは他のプラットフォームでコンパイルして動作します。



まず、サウンドファイルを生成するプログラムを作成します。 wavファイルは最も単純なものとして扱われました。 ここでその構造について読んでください

要するに、wavファイルの構造は次のように説明されています:ファイル形式を記述するヘッダー、そして(この場合)長さの16ビットデータ(ポイントされた)の配列があります:Sample_frequency * t秒または44100 * t個。



オーディオファイルを作成するために、 ここで例を取り上げました 。 私はそれを少し修正し、エラーを修正し、編集した最終バージョンはここのgithubにあります



github.com/dlinyj/generate_wav



100 Hzの純粋な正弦周波数を持つ2秒のサウンドファイルを生成します。 これを行うには、次の方法でプログラムを変更します。



 #define S_RATE (44100) //  #define BUF_SIZE (S_RATE*10) /* 2 second buffer */ …. int main(int argc, char * argv[]) { ... float amplitude = 32000; //    float freq_Hz = 100; //  /* fill buffer with a sine wave */ for (i=0; i<BUF_SIZE; i++) { buffer[i] +=(int)(amplitude * sin((float)(2*M_PI*i*freq_Hz/S_RATE))); } write_wav("test.wav", BUF_SIZE, buffer, S_RATE); return 0; }
      
      







純粋な正弦波の式は上記で説明した式に対応するという事実に注意を促します。 32000の振幅(32767を取ることができます)は、16ビット数(マイナス32767からプラス32767まで)を取ることができる値に対応します。



その結果、次のファイルが得られます (任意の音声再生プログラムで聞くこともできます )。 このaudacityファイルを開き、信号グラフが実際に純粋な正弦に対応することを確認します。





純粋なチューブサイン



このサインのスペクトルを見てみましょう(分析->スペクトルグラフの作成)





スペクトルグラフ



100 Hz(対数目盛)で明確なピークが見られます。 スペクトルとは何ですか? これは振幅周波数応答です。 位相周波数特性もあります。 覚えている場合は、信号を作成するには、その周波数、振幅、および位相を知る必要があると上記で述べました。 したがって、信号からこれらのパラメーターを取得できます。 この場合、周波数と振幅の対応関係のグラフがあり、振幅は実単位ではなくデシベル単位です。







デシベルで表された値は、同じ名前の物理量に対する物理量の無次元比の10進数の対数に数値的に等しく、これは最初の1に10を掛けたものです。



この場合、振幅の対数に10を掛けただけです。対数目盛は、信号を処理するときに使用すると便利です。



正直に言うと、私はこのプログラムのスペクトルアナライザーがあまり好きではないので、特に難しいことではないので、 ブラックジャックと売春婦で自分書くことにしました。



スペクトルアナライザーを作成します





ここは退屈かもしれませんので、次の章に進んでください。



私はコードの足跡を投稿する意味がないことを完全によく理解しているので、本当に興味がある人は自分で見つけて掘り下げ、興味がない人は退屈するでしょう、私はwavファイルスペクトルアナライザを書くことの主なポイントにのみ焦点を当てます。



最初に、wavファイルを読み取る必要があります。 そこで、ヘッダーを読んで、このファイルに含まれるものを理解する必要があります。 このファイルを読み込むためのオプションの海を実装しませんでしたが、1つだけで解決しました。 ここからファイルを読み取る例はほとんど変更なしで取得されましたが、私見は素晴らしい例です。 Python実装もあります。



次に必要なのは、高速フーリエ変換です。 これは、点の有限セットから元の信号のベクトルaを取得できるようにするまさに変換です。 まだあなたを怖がらせないでください、私はさらに説明します。

繰り返しますが、彼は自転車を発明しませんでしたが、 ここから既製の例を取りました。



プログラムの動作を説明するためには、高速フーリエ変換とは何かを説明する必要があることを理解しています。これは少なくとも1つの非酸性の記事です。



まず、配列を分離します。



  c = calloc(size_array*2, sizeof(float)); //    in = calloc(size_array*2, sizeof(float)); //  out = calloc(size_array*2, sizeof(float)); // 
      
      







プログラムでは、size_arrayの長さの配列(wavファイルのヘッダーから取得)にデータを読み取ります。



  while( fread(&value,sizeof(value),1,wav) ) { in[j]=(float)value; j+=2; if (j > 2*size_array) break; }
      
      







高速フーリエ変換の配列は、シーケンス{re [0]、im [0]、re [1]、im [1]、... re [fft_size-1]、im [fft_size-1]}である必要があります。ここで、fft_size = 1 << pは、FFTポイントの数です。 私は通常の言語で説明します:

複素数の配列です。 複素フーリエ変換がどこで使用されるか想像するのも怖いですが、この場合、虚数部はゼロで、実数部は配列の各ポイントの値に等しくなります。

高速フーリエ変換のもう1つの特徴は、2のべき乗のみの倍数である配列を計算することです。 その結果、2の最小電力を計算する必要があります。



 int p2=(int)(log2(header.bytes_in_data/header.bytes_by_capture));
      
      







データのバイト数を1ポイントのバイト数で割った対数。



その後、ターニングファクターを考慮します。



 fft_make(p2,c);//       (   ,     ).
      
      







そして、読み取り配列をフーリエ変換器に送ります。



 fft_calc(p2, c, in, out, 1); //( ,     ).
      
      







出力では、{re [0]、im [0]、re [1]、im [1]、... re [fft_size-1]、im [fft_size-1]}の形式の複素数を取得します。 複素数が何であるかを知らない人のために、私は説明します。 たくさんの回転ベクトルとたくさんのgifを使ってこの記事を始めたのは、何の理由もありませんでした。 したがって、複素平面上のベクトルは、実座標a1と虚座標a2によって決まります。 または、長さ(これは振幅Am)と角度Psi(位相)です。





複素平面ベクトル



size_array = 2 ^ p2に注意してください。 配列の最初のポイントは0 Hz(一定)の周波数に対応し、最後のポイントはサンプリング周波数、つまり44100 Hzに対応します。 その結果、各ポイントに対応する周波数を計算する必要がありますが、これは周波数デルタによって異なります。



 double delta=((float)header.frequency)/(float)size_array; //    .
      
      







振幅の配列の配置:



  double * ampl; ampl = calloc(size_array*2, sizeof(double));
      
      







そして画像を見てください:振幅はベクトルの長さです。 そして、実軸と虚軸に投影しています。 その結果、直角三角形になり、ここでピタゴラスの定理を思い出し、各ベクトルの長さを考慮して、すぐにテキストファイルに書き込みます。



 for(i=0;i<(size_array);i+=2) { fprintf(logfile,"%.6f %f\n",cur_freq, (sqrt(out[i]*out[i]+out[i+1]*out[i+1]))); cur_freq+=delta; }
      
      





その結果、この種類のファイルを取得します。



 … 11.439514 10.943008 11.607742 56.649738 11.775970 15.652428 11.944199 21.872342 12.112427 30.635371 12.280655 30.329171 12.448883 11.932371 12.617111 20.777617 ...
      
      







プログラムの最終バージョンは、次のgithubにあります。

github.com/dlinyj/fft



やってみます!





次に、サインファイルを鳴らす結果のプログラムをフィードします



 ./fft_an ../generate_wav/sin\ 100\ Hz.wav format: 16 bits, PCM uncompressed, channel 1, freq 44100, 88200 bytes per sec, 2 bytes by capture, 2 bits per sample, 882000 bytes in data chunk=441000 log2=18 size array=262144 wav format Max Freq = 99.928 , amp =7216.136
      
      







そして、周波数応答のテキストファイルを取得します。 sluplotの助けを借りてチャートを作成する



構築するスクリプト:



 #! /usr/bin/gnuplot -persist set terminal postscript eps enhanced color solid set output "result.ps" #set terminal png size 800, 600 #set output "result.png" set grid xtics ytics set log xy set xlabel "Freq, Hz" set ylabel "Amp, dB" set xrange [1:22050] #set yrange [0.00001:100000] plot "test.txt" using 1:2 title "AFC" with lines linestyle 1
      
      







X:set xrange [1:22050]のポイント数に関するスクリプトの制限に注意してください。 44100のサンプリング周波数があり、コテルニコフの定理を思い出せば、信号周波数はサンプリング周波数の半分より高くすることはできません。したがって、22050 Hzを超える信号は対象になりません。 なぜそうなのか、専門文献を読むことをお勧めします。

だから(ドラムロール)、スクリプトを実行して確認してください:





信号のスペクトル



100 Hzの鋭いピークに注意してください。 軸が対数目盛であることを忘れないでください! 私が思うに、右側の羊毛はフーリエ変換誤差です(ここでウィンドウが記憶されます)。



そして、甘やかしましょうか?





さあ! 他の信号のスペクトルを見てみましょう!



ノイズの周りに...


まず、ノイズスペクトルを作成します。 トピックは、ノイズ、ランダム信号などです。 別のコースに値する。 しかし、軽く触れます。 wavファイルを生成するためにプログラムを変更し、1つの手順を追加します。



 double d_random(double min, double max) { return min + (max - min) / RAND_MAX * rand(); }
      
      







指定された範囲で乱数を生成します。 その結果、メインは次のようになります。



 int main(int argc, char * argv[]) { int i; float amplitude = 32000; srand((unsigned int)time(0)); //    for (i=0; i<BUF_SIZE; i++) { buffer[i] +=(int)amplitude*d_random(-1.0, 1.0); //nois } write_wav("test.wav", BUF_SIZE, buffer, S_RATE); return 0; }
      
      







ファイルを生成します (リスニングをお勧めします )。 それを大胆に見てみましょう。





大胆なシグナル



audacityプログラムのスペクトルを見てみましょう。





スペクトル



そして、プログラムを使用してスペクトルを見てみましょう。





私たちの範囲



私は非常に興味深い事実とノイズの特徴に注目したいと思います-それはすべての高調波のスペクトルを含んでいます。 グラフからわかるように、スペクトルは非常に均一です。 原則として、 ホワイトノイズは帯域幅の周波数分析、たとえばオーディオ機器に使用されます。 他の種類のノイズがあります: ピンク、ブルー、その他 。 宿題-それらがどのように異なるかを調べます。



コンポートはどうですか?




次に、別の興味深い信号、蛇行を見てみましょう。 上記のフーリエ級数のさまざまな信号の拡大の兆候を引用しましたが、蛇行がどのようにレイアウトされているかを確認し、それを紙に書いて、続けます。



周波数25 Hzの方形波を生成するために、wavファイルジェネレーターを再度変更します。



 int main(int argc, char * argv[]) { int i; short int meandr_value=32767; /* fill buffer with a sine wave */ for (i=0; i<BUF_SIZE; i++) { //meandr if (!(i%(S_RATE/((int)freq_Hz/2)))) { if (meandr_value==32767) { meandr_value=-32767; } else { meandr_value=32767; } } buffer[i]=meandr_value; } write_wav("test.wav", BUF_SIZE, buffer, S_RATE); return 0; }
      
      







その結果、サウンドファイルを取得します (もう一度、聞くことをお勧めします)。すぐにaudacityで見る必要があります。





His下-健康な人の蛇行または蛇行



苦労せずに、そのスペクトルを見てみましょう。





スペクトル蛇行



これまでのところ、何かがあまり明確ではありません...そして、最初のいくつかの高調波を見てみましょう:





第一高調波



完全に異なる問題! さて、タブレットを見てみましょう。 見て、1、3、5などしかありません。 奇数倍音。 25 Hzの最初の高調波、次の(3番目の)75 Hz、125 Hzなどがありますが、振幅は徐々に減少しています。 理論は実践に収束します!

今注目! 実際には、私たちの蛇行信号にはますます高い周波数の無限の高調波がありますが、原則として、実際の電気回路は特定の周波数を超える周波数を通過させることはできません(トラックのインダクタンスとキャパシタンスのため)。 その結果、オシロスコープ画面に次の信号が表示されることがよくあります。





喫煙者の蛇行



この写真は、 ウィキペディアの写真に似ています。蛇行の例では、すべての周波数ではなく、最初の数個のみが撮影されています。





最初の高調波の合計と信号の変化



蛇行は無線工学でも積極的に使用されています(これはすべてのデジタルテクノロジーの基礎であると言わなければなりません)。長い回路では、母親が認識しないようにフィルターで除去できることを理解する価値があります。 また、さまざまなデバイスの周波数応答を確認するためにも使用されます。 別の興味深い事実は、超小型回路自体が数十MHzの蛇行を生成し、その高次高調波がテレビの周波数で数百MHzの周波数を持ち、高次高調波がテレビの放送信号を正常に減衰させるときに、TVジャマーが高次高調波の原理で動作したことです。



一般に、そのような実験のテーマは無限であり、あなたは今、あなた自身でそれを続けることができます。



読書の推奨事項











私たちがしていることを理解していない人、またはその逆、理解しているがさらに理解したい人、およびDSPを勉強している学生のために、この本を強くお勧めします。 これは、この投稿の著者であるダミー用DSPです。 そこでは、子供にとってさえアクセスしやすい言語が最も複雑な概念を話します。



おわりに





結論として、数学は科学の女王ですが、実際の応用がなければ、多くの人が数学に興味を失います。 この投稿が、信号処理や一般的なアナログ回路(脳が流れ出ないように耳をふさぐ!)などのすばらしいテーマを学ぶことを奨励することを願っています。 :)

頑張って



All Articles