前回は、ブログに数式を記入するだけで力尽き、本文まで書けなかったので、今回はその続きです。
実際にC#でプログラムしてみました。 普通は、フーリエ変換というとFFTだと思うのですが、データ数が少なく、リアルタイムで処理するわけでもないので、式通りの離散的フーリエ変換を行いました。 こっちの方が断然簡単です。
いきなり、FPGAの信号を処理してもいいのですが、うまくいかなかった時にデバッグを簡単にするために、とりあえずダミーのサイン波を入力して、表示させてみました。
DFTの結果は、複素数になるので中央より左側が負の周波数領域になります。中央が40MHzで、左端、右端とも0MHzになります。 少し見辛いですが、5MHzのところに赤い線で表示した信号が出ています。
ソースです。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 |
public partial class Form1 : Form { bool sindata_ready = false; bool sadata_ready = false; double[] sin_data = new double[1024]; double[] mag = new double[1024]; Complex[] Xk; Speana dg1; DrawGraph dg2; public Form1() { InitializeComponent(); } private void button1_Click(object sender, EventArgs e) { double sig_freq = Double.Parse(textBox1.Text); double samplerate = Double.Parse(textBox2.Text); double magnitude = Double.Parse(textBox3.Text); //テスト用sin波作成 for(int i = 0; i < 1024; i++ ) { sin_data[i] = magnitude * Math.Sin(i * 2.0 * Math.PI * sig_freq / samplerate); } sindata_ready = true; Xk = new Complex[sin_data.Length]; //Dft実行 Xk = DftExecute(sin_data); //dB変換 1e-9は、0になるのを避けるため for (int i = 0; i < Xk.Length; i++) mag[i] = Math.Log10(Xk[i].Magnitude+1e-9) * 20; sadata_ready = true; pictureBox1.Refresh(); pictureBox2.Refresh(); } private void pictureBox1_Paint(object sender, PaintEventArgs e) { dg1.DrawAxis(e.Graphics); if (sadata_ready) dg1.DrawSpectrum(e.Graphics, mag); } private void pictureBox2_Paint(object sender, PaintEventArgs e) { dg2.DrawAxis(e.Graphics); if (sindata_ready) dg2.DrawData(e.Graphics, sin_data); } private void Form1_Load(object sender, EventArgs e) { dg1 = new Speana(pictureBox1); dg2 = new DrawGraph(pictureBox2); } Complex[] DftExecute(double[] xn) { Complex[] Xk = new Complex[xn.Length]; for (int i = 0; i < xn.Length; i++) { Xk[i] = 0; for (int j = 0; j < xn.Length; j++) { Xk[i] = Xk[i] + xn[j] * Complex.Exp(2 * Math.PI * j * i * Complex.ImaginaryOne / xn.Length); } } return Xk; } } |
ところどころにコメントをいれましたが、最後のDftExecuteがフーリエ変換の部分です。 前回記載した数式の通りに計算しています。
波形表示のクラスです。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 |
class DrawGraph { private PictureBox picturebox_; private float x0_; //x軸原点 private float y0_; //y軸原点 private float y1_; //y軸1に対するピクセル値 private float width_; private Font myfont_ = new Font("MS UI Gothic", 10); private SolidBrush mysolidbrush_ = new SolidBrush(Color.Black); private const float x_offset = 10.0F; private const float y_offset = 10.0F; //コンストラクタ public DrawGraph(PictureBox myPictureBox) { picturebox_ = myPictureBox; x0_ = picturebox_.Left; y0_ = picturebox_.ClientSize.Height - 1; y1_ = y0_ /(float)(1.2* 1024); width_ = picturebox_.ClientSize.Width - 1; } public void DrawAxis(Graphics gr) { DrawVerticalLines(gr, Pens.Black); DrawHolizontalLines(gr, Pens.Black); } private void DrawHolizontalLines(Graphics gr, Pen mypen) { float yn = y0_ / 10; for (int i = 0; i <= y0_; i += 1) gr.DrawLine(mypen, 0, y0_ - yn * i, width_, y0_ - yn * i); } private void DrawVerticalLines(Graphics gr, Pen mypen) { float xn = width_ / 10; for (int i = 0; i <= (int)width_; i += 1) gr.DrawLine(mypen, xn * i, 0, xn * i, y0_); } public void DrawData(Graphics gr, double[] gdata) { Pen mypen = new Pen(Color.Red); PointF oldPoint = new PointF(0, (float)(y0_ - y1_ * gdata[0] - y0_/2)); float tick = width_ / gdata.Length*5; for (int i = 1; i < gdata.Length - 1; i++) { mypen.Color = Color.Red; PointF newPoint = new PointF(i * tick, (float)(y0_ - y1_ * gdata[i] - y0_/2)); gr.DrawLine(mypen, oldPoint, newPoint); oldPoint = newPoint; } } } |
スペクトル表示用のクラスです。 縦軸はdB表示にしていますが、入力信号がダミーなので適当です。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 |
class Speana { private PictureBox picturebox_; private float x0_; //x軸原点 private float y0_; //y軸原点 private float width_; private const float x_offset = 10.0F; private const float y_offset = 10.0F; public Speana(PictureBox myPitureBox) { picturebox_ = myPitureBox; x0_ = picturebox_.Left; y0_ = picturebox_.ClientSize.Height - 1; width_ = picturebox_.ClientSize.Width - 1; } public void DrawAxis(Graphics gr) { float step = width_ / 16; for (int i = 0; i <= 16; i++) gr.DrawLine(Pens.Black, step * i, 0, step * i, y0_); step = y0_ / 10; for (int i = 0; i <= 10; i++) gr.DrawLine(Pens.Black, 0, y0_ - step * i, width_, y0_ - step * i); step = step/2; for (int i = 0; i <= 20; i++) gr.DrawLine(Pens.Black, width_/2, y0_ - step * i, width_/2 + 3, y0_ - step * i); } public void DrawSpectrum(Graphics gr, double[] data) { float xtick = width_ / 1024f; PointF oldPoint = new PointF(0, y0_- (float)(data[0])); PointF newPoint = new PointF(); for (int i = 1; i < data.Length; i++) { newPoint.X = i * xtick; newPoint.Y = y0_-(float)data[i]; gr.DrawLine(Pens.Red, oldPoint, newPoint); oldPoint = newPoint; } } } |
DFTのプログラムが正常に動作していることが確認できたので、次回は実際のFPGAから入力した信号でやってみたいと思います。