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周波数解析のプログラムテストのため、2個のデーターを簡単に出力できるようにしました。
一周期の分割数、波形の数の入力ボックスの下の矢印で、値を設定すると、2個のデーター数になります。
矢印を使用せずに、自由な値を設定することもできます。
ファイル作成にあるサンプリング時間は、周波数解析用として、単にファイルに書き込む為のものです。

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.

    download Waveform.zip

各種プログラム計算例に戻る

      最初に戻る