2020/05/01 フーリエ変換による周波数解析のテストに使用できるように、ファイル出力を追加しました。
波形と数式
振動波形を、プログラム上で再現する為の計算式です。
強制振動の加速度に使用すれば、色々な強制波形に対する応答のグラフを得ることが出来ます。
周波数解析の為のサンプルデーターの作成も出来ます。
波形の分割数が波形の一サイクルのデーターの数ですが、分割数を少なくすると波形の再現が出来なくなり、直線や三角波になります。、
正弦波形
Sin(ωt)
此処ではCos(ωt+α) α=pi/2 としてSin(ωt)と同じ計算をしています。
ノコギリ波
ノコギリ波には、方向があります。
①は、傾斜が右下がりで、②は周期がαのノコギリ波。
③は、フーリエ級数によるノコギリ波です。
右上がりのノコギリ波にする場合は、符号を反転するだけです。
矩形波
①②矩形波の数式で、a,b は high と low の比率となります。
③はフーリェ級数による矩形波です。
三角波
①階段波
③矩形波
④ノコギリ波
⑤は、三角波を表す数式で
①から④の数式を利用しています。
a
の値は0.5にとります。
⑥が一番簡単な三角波となる数式です。
⑦は、フーリェ級数による三角波です。
階段波
①②③は、階段波を表す数式です。
階段波を表すフーリエ級数は複雑になるのでありません。
ファイルへの出力
8192 1.00000000000000E+0000 <----------------
データー数 サンプリング時間
-5.00000000000000E-0001 <---------------- データー
-4.92187500000000E-0001
-4.84375000000000E-0001
-------------
-------------
-------------
ファイルへの出力は、最初に、データー数とサンプリング時間を書き込み、
次にデーター数分の値を書き込みます。
フーリエ変換のプログラムで読み込んで使用できます。
プログラム
一般数式による波形と、フーリェ級数による波形の出力はラジオボタンで選択します。
フーリェ級数による場合は、1サイクルを表示する為の周波数の数があり、周波数の数が多いほどきれいな波形となります。
FFT周波数解析のプログラムテストのため、2n個のデーターを簡単に出力できるようにしました。
一周期の分割数、波形の数の入力ボックスの下の矢印で、値を設定すると、2n個のデーター数になります。
矢印を使用せずに、自由な値を設定することもできます。
ファイル作成にあるサンプリング時間は、周波数解析用として、単にファイルに書き込む為のものです。
unit Main; interface uses Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics, Vcl.Controls, Vcl.Forms, Vcl.Dialogs, VclTee.TeeGDIPlus, VCLTee.TeEngine, VCLTee.Series, Vcl.ExtCtrls, VCLTee.TeeProcs, VCLTee.Chart, Vcl.StdCtrls, Vcl.Samples.Spin, Vcl.ComCtrls ,System.UITypes; type TWaveForm1 = class(TForm) GroupBox1: TGroupBox; RadioButton1: TRadioButton; RadioButton2: TRadioButton; Chart1: TChart; Series1: TFastLineSeries; RadioButton3: TRadioButton; aEdit: TLabeledEdit; bEdit: TLabeledEdit; RadioButton4: TRadioButton; RadioButton5: TRadioButton; GroupBox2: TGroupBox; RadioButton6: TRadioButton; MEdit: TLabeledEdit; RadioButton7: TRadioButton; RadioButton8: TRadioButton; RadioButton9: TRadioButton; GroupBox3: TGroupBox; WaveButton: TRadioButton; anyButton: TRadioButton; ndEdit: TLabeledEdit; TEdit: TLabeledEdit; UpDown1: TUpDown; UpDown2: TUpDown; MakeBtn: TButton; GroupBox4: TGroupBox; smpTEdit: TLabeledEdit; SaveBtn: TButton; ReadBtn: TButton; SaveDialog1: TSaveDialog; OpenDialog1: TOpenDialog; OffsetEdit: TLabeledEdit; RadioButton10: TRadioButton; procedure FormCreate(Sender: TObject); procedure UpDown1MouseUp(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer); procedure UpDown2MouseUp(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer); procedure MakeBtnClick(Sender: TObject); procedure SaveBtnClick(Sender: TObject); procedure ReadBtnClick(Sender: TObject); private { Private 宣言 } procedure Wavemake; procedure Any_waveMake(M, L: integer); procedure anyWave; function DataSave: boolean; public { Public 宣言 } end; var WaveForm1: TWaveForm1; implementation {$R *.dfm} uses math; var N : integer; // データー数 TS : double; // サンプリング時間 nd : integer; // 一周期の分割数 T : integer; // 繰り返し数 FS : TextFile; // テキストファイル SF : boolean; // ファイルSaveフラグ offset : double; // オフセット //============================== 数式 ============================ // 階段波形 function f1(x: double): double; begin result := floor(x); end; // 正弦波 function c1(x: double): double; begin result := cos(2 * pi * x - pi / 2); end; // のこぎり波形 function s1(x: double): double; begin result := x - floor(x); end; //---------------------------------------------------- // 階段波 function fl(x: double): double; begin result := (floor(x) - floor(-x) - 1) / 2; end; function k(a, b, x: double): double; begin result := fl((x + a)/(a + b)) - fl(x / (a + b)); end; // 矩形波 kk(a, b, x: double); function KK(a, b, x: double): double; begin result := 2 * k(a,b,x) - 1; end; // ノコギリ波 function s2(x: double): double; begin result := x + floor(-x) + 1; end; // 三角波 function s(a, x: double): double; begin result := 2 * abs(s2(x / a) - 1 / 2) * kk(a, a, x + a / 2); end; //------------------------------------------------------------ procedure TWaveForm1.Wavemake; var dt : double; x : double; y : double; k : integer; a, b: double; ymax, ymin : double; begin y := 0; a := 0.5; b := 1; if RadioButton3.Checked then begin val(aEdit.text, a, k); if (k <> 0) or (a <= 0) then begin application.MessageBox('aの値に間違いがあります。','注意',0); exit; end; val(bEdit.text, b, k); if (k <> 0) or (b <= 0) then begin application.MessageBox('bの値に間違いがあります。','注意',0); exit; end; end; val(ndEdit.text, nd, k); if (k <> 0) or (nd < 1) then begin application.MessageBox('1周期の分割数の値に間違いがあります。','注意',0); exit; end; N := T * nd; // データーの数 dt := 1 / nd; // 1波形の分割数 // n := round(t / dt); // データーの数 if SF then Writeln(FS, N, TS); ymax := 0; ymin := 0; Series1.Clear; for k := 0 to N do begin // 作図 x := dt * k; if RadioButton1.Checked then y := c1(x); // 正弦波 if RadioButton2.Checked then y := s1(x) * 2 - 1; // 逆ノコギリ波 if RadioButton3.Checked then y := KK(a, b, x); // 矩形波 // if RadioButton3.Checked then y := (f1((x + a)/(a + b) - f1(x/(a + b)))) * 2 - 1; // 矩形波 if RadioButton4.Checked then y := (1 - s1(x)) * 2 - 1; // ノコギリ波 if RadioButton5.Checked then y := s(0.5, x); // 三角波 // if RadioButton5.Checked then y := (arccos(cos(x * 2 * pi))) / pi * 2 - 1; // 三角波 if K < N then if RadioButton10.Checked then y := f1(x); // 階段波 y := y + offset; Series1.AddXY(x, y); if SF and (K < N) then Writeln(FS, y); if y < ymin then ymin := y; if y > ymax then ymax := y; end; Chart1.LeftAxis.Maximum := ymax + 1; Chart1.LeftAxis.Minimum := ymin - 1; end; //==================== フーリエ級数 ============================================ // 矩形波 function Square_wave(t: double; K: integer): double; var omega : double; A, I : integer; S : double; begin omega := 2 * pi; s := 0; A := 1; for I := 1 to K do begin S := S + sin(A * omega * t) / A; A := A + 2; end; result := S * 4 / pi; end; // のこぎり波 function sawtooth_wave(t: double; K: integer): double; var omega : double; A, I : integer; S : double; begin omega := 2 * pi; s := 0; A := 1; for I := 1 to K do begin S := S + sin(A * omega * t) / A; inc(A); end; result := S * 2 / pi; end; // 三角波 function triangle_wave(t: double; K: integer): double; var omega : double; A, I : integer; S : double; begin omega := 2 * pi; s := 0; A := 1; for I := 1 to K do begin S := S + sin(A * omega * t) / (A * A) * sin(A * pi / 2); inc(A); end; result := S * 8 / (pi * pi); end; // フーリエ級数による波形 procedure TWaveForm1.Any_waveMake(M, L: integer); var dt, ft :double; k : integer; y : double; ymax, ymin : double; begin y := 0; N := L * T; dt := 1 / L; if SF then Writeln(FS, N, TS); ymax := 0; ymin := 0; Series1.Clear; for k := 0 to N do begin ft := k * dt; if RadioButton6.Checked then y := Square_wave(ft, M); if RadioButton7.Checked then y := sawtooth_wave(ft, M); if RadioButton8.Checked then y := -sawtooth_wave(ft, M); if RadioButton9.Checked then y := triangle_wave(ft, M); y := y + offset; Series1.AddXY(ft, y); if SF and (K < N) then Writeln(FS, y); if y < ymin then ymin := y; if y > ymax then ymax := y; end; Chart1.LeftAxis.Maximum := ymax + 1; Chart1.LeftAxis.Minimum := ymin - 1; end; procedure TWaveForm1.anyWave; var M, C: integer; begin val(MEdit.text, M, C); if (C <> 0) or (M < 1) then begin application.MessageBox('正弦波の数。','注意',0); exit; end; Any_waveMake(M, nd); end; var ndi, Ti : integer; //================== Data 作成 ======================= procedure TWaveForm1.MakeBtnClick(Sender: TObject); var ch : integer; begin val(TEdit.Text, T, ch); if ch <> 0 then begin application.MessageBox('波形の数に間違いがあります。','',0); exit; end; if T < 1 then begin application.MessageBox('波形の数が少なすぎます。','',0); exit; end; val(ndEdit.Text, nd, ch); if ch <> 0 then begin application.MessageBox('分割数に間違いがあります。','',0); exit; end; if nd < 2 then begin application.MessageBox('一周期の分割数が少なすぎます','',0); exit; end; if (anyButton.Checked) and (nd < 4) then begin application.MessageBox('Fourierの場合一周期の分割数が少なすぎます','',0); exit; end; val(offsetEdit.Text, offset, ch); if ch <> 0 then begin application.MessageBox('オフセットに間違いがあります。','注意', 0); exit; end; if WaveButton.Checked then Wavemake; if anyButton.Checked then anyWave; end; //----------------------- 分割数設定 ----------------------------- procedure TWaveForm1.UpDown1MouseUp(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer); var I : integer; begin if UpDown1.Position > 8 then UpDown1.Position := 8; if UpDown1.Position < 1 then UpDown1.Position := 1; ndi := 1; for I := 1 to UpDown1.Position do ndi := ndi * 2; ndEdit.Text := intTostr(ndi); nd := ndi; end; //------------------------ 繰り返し数設定 ------------------------------ procedure TWaveForm1.UpDown2MouseUp(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer); var I : integer; begin if UpDown2.Position > 6 then UpDown2.Position := 6; if UpDown2.Position < 1 then UpDown2.Position := 1; Ti := 1; for I := 1 to UpDown2.Position do Ti := Ti * 2; TEdit.Text := intTostr(Ti); T := Ti; end; //======================== ファイル操作 ================================ function TWaveForm1.DataSave: boolean; var ch : integer; begin result := False; val(smpTEdit.Text, TS, ch); if ch <> 0 then begin application.MessageBox('サンプリングの値に間違いがあります。','注意',0); exit; end; val(offsetEdit.Text, offset, ch); if ch <> 0 then begin application.MessageBox('オフセットに間違いがあります。','注意', 0); exit; end; if (anyButton.Checked) and (nd < 4) then begin application.MessageBox('Fourierの場合一周期の分割数が少なすぎます','',0); exit; end; result := True; end; procedure TWaveForm1.SaveBtnClick(Sender: TObject); begin SF := False; if not DataSave then exit; SaveDialog1.Options := [ofHideReadOnly,ofNoReadOnlyReturn]; ReadBtn.Enabled := False; SaveBtn.Enabled := False; if SaveDialog1.Execute then begin if AnsiLowerCase(ExtractFileExt(SaveDialog1.FileName)) <> '.fre' then SaveDialog1.FileName := SaveDialog1.FileName + '.fre'; if FileExists(SaveDialog1.FileName) then begin if MessageDlg('既にファイルが存在します 上書きしますか',mtConfirmation, [mbYes, mbNo], 0) = mrNo then exit; end; try AssignFile(FS, SaveDialog1.FileName); ReWrite(FS); SF := True; MakeBtnClick(nil); finally CloseFile(FS); SF := False; end; end; ReadBtn.Enabled := True; SaveBtn.Enabled := True; end; procedure TWaveForm1.ReadBtnClick(Sender: TObject); var KK : integer; inStr : string; data : double; fdata : double; begin ReadBtn.Enabled := False; SaveBtn.Enabled := False; if OpenDialog1.Execute then begin Series1.Clear; try AssignFile(FS,OpenDialog1.FileName); Reset(FS); Readln(FS, N, TS); smpTEdit.Text := floatTostr(TS); Readln(FS, inStr); fdata := strTofloat(inStr); Series1.AddXY(0, fdata); for KK := 1 to N -1 do begin Readln(FS, inStr); data := strTofloat(inStr); Series1.AddXY(KK, data); end; Series1.AddXY(N, fdata); finally CloseFile(FS); end; end; ReadBtn.Enabled := True; SaveBtn.Enabled := True; end; //========================= 初期設定 ========================= procedure TWaveForm1.FormCreate(Sender: TObject); var I : integer; begin UpDown1.Position := 7; ndi := 1; for I := 1 to UpDown1.Position do ndi := ndi * 2; ndEdit.Text := intTostr(ndi); nd := ndi; Ti := 1; UpDown2.Position := 2; for I := 1 to UpDown2.Position do Ti := Ti * 2; T := Ti; SF := False; end; end.