// 重み付きチェビシェフ近似による直線位相 FIR フィルタの設計
// （Parks-McClellan のアルゴリズムによる）
//      T. W. Parks and C. S. Burrus: “Digital filter design,”
//      Chapter 3, Joh Wiley & Sons, Inc., 1987. 参照
// 設計結果の表示スペースをマウスにより拡大可能にしたバージョン
//
// Embarcadero C++ Builder 2010 で再ビルド  2010/02/23
//---------------------------------------------------------------------------
#include <vcl.h>
#include <cmath>
using namespace std;

#include "Unit1.h"
#include "Unit2.h"
#include "MK_misc.hpp"
#pragma hdrstop

#include "Remez.hpp"

const int fDiv = 1024;                  // 振幅特性描画の際の点数
const double Amin = 1.0e-10;            // dB 値の最小値に対応する値
const double AmindB = 20.0*log10(Amin); // dB 値の最小値

VectorD hk(MAX_TAPS);                   // 係数
VectorD LogA_k(fDiv+1);                 // 振幅特性の値
VectorD bands(2*MAX_BANDS);             // 帯域端周波数

void inline Line(TCanvas *Cv, int x1, int y1, int x2, int y2)
{ Cv->MoveTo(x1, y1); Cv->LineTo(x2, y2); }

//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
    : TForm(Owner)
{
    const int x0 = 176;
    const int y0 = Form1->ClientHeight - 97;

    X0d =  50;          // x軸の原点（周波数 = 0）
    Y0d =  26;          // y軸の原点（0 dB）
    Xfr = 105;          // x軸の1目盛り（(fs/2)/5）
    YdB =  60;          // y軸の1目盛り（表示dBの幅/8）

    SetEdit(pEditWght, MAX_BANDS, x0, y0+24*3, "1");
    SetEdit(pEditGain, MAX_BANDS, x0, y0+24*2, "1");
    for (int i=1; i<MAX_BANDS; i++) pEditGain[i]->Text = "0";
    SetEdit(pEditBdH,  MAX_BANDS, x0, y0+24, FloatToStr(StrToFloat(Edit2->Text)/2));
    SetEdit(pEditBdL,  MAX_BANDS, x0, y0,    "0");
    pEditBdH[0]->Text = "1.0";
    pEditBdL[1]->Text = "1.5";
    for( int i=0; i<MAX_BANDS; i++ )
    {
        pLabelBnd[i] = new TLabel(this);
        pLabelBnd[i]->Parent = this;
        pLabelBnd[i]->Width = 10;
        pLabelBnd[i]->Height = 12;
        pLabelBnd[i]->Top = y0 - 20;
        pLabelBnd[i]->Left = x0 + i*88 + 18;
        pLabelBnd[i]->Caption = (AnsiString)"帯域 "+IntToStr(i+1);
        if (i<2) pLabelBnd[i]->Enabled = true;
        else     pLabelBnd[i]->Enabled = false;
    }
    Label3->Top = y0+3;
    Label4->Top = Label3->Top + 24;
    Label5->Top = Label4->Top + 24;
    Label6->Top = Label5->Top + 24;
    UpDown1->Max = MAX_BANDS;
    ComboBox1->ItemIndex = 4;                   // 表示範囲を -80 dB までに設定

    RichEdit1->Top = y0 - 386;
    width_of_RichEdit1  = RichEdit1->Width;     // 起動時の RichEdit1 の幅を保存
    height_of_RichEdit1 = RichEdit1->Height;    // 起動時の RichEdit1 の高さを保存
    width_of_Form1  = Form1->Width;             // 設計時の Form1 の幅
    height_of_Form1 = Form1->Height;            // 設計時の Form1 の高さ

    calc_ok = false;
}
//---------------------------------------------------------------------------

__fastcall TForm1::~TForm1()
{
    for( int i=0; i < MAX_BANDS; i++ ) delete pEditBdL[i];
    for( int i=0; i < MAX_BANDS; i++ ) delete pEditBdH[i];
    for( int i=0; i < MAX_BANDS; i++ ) delete pEditGain[i];
    for( int i=0; i < MAX_BANDS; i++ ) delete pEditWght[i];
    for( int i=0; i < MAX_BANDS; i++ ) delete pLabelBnd[i];
}
//---------------------------------------------------------------------------

// “設計”ボタン
void __fastcall TForm1::Button1Click(TObject *Sender)
{
    AnsiString buf, buf2;
    VectorD weights(MAX_BANDS), gain(MAX_BANDS);
    VectorD deviat(MAX_BANDS), dev_dB(MAX_BANDS);

    if (!InputSpec(weights, gain)) return;      // 設計パラメータの入力
    Button1->Caption = "計算中";
    Button1->Update();
    RichEdit1->Lines->Clear();
    if (!Remez(nTaps, nBands, bands, gain, weights, deviat, jtype, fs))
    {
        ShowMessage("収束しません．パラメータの与え方を変えてください．");
        return;
    }

// 係数の出力
    for (int j=0; j<nTaps/4; j++)
    {
        buf.sprintf("        %15.8e, %15.8e, %15.8e, %15.8e,",
                         hk[4*j], hk[4*j+1], hk[4*j+2], hk[4*j+3]);
        if ( ((nTaps % 4) == 0) && (j == nTaps/4-1) )
        {
            buf.Delete(buf.Length(), 1);
            buf = buf + " };";
        }
        RichEdit1->Lines->Add(buf);
    }
    if ((nTaps % 4) != 0)
    {
        buf = "        ";
        for (int j=(nTaps/4)*4; j<nTaps; j++)
            buf = buf + buf2.sprintf("%15.8e, ", hk[j]);
        buf.Delete(buf.Length()-1, 2);
        buf = buf + " };";
        RichEdit1->Lines->Add(buf);
    }
// 設計パラメータ及びその他の結果の出力
    RichEdit1->Lines->Add("//");
    RichEdit1->Lines->Add("// 次数：　　　　 "+ Edit1->Text);
    RichEdit1->Lines->Add("// 標本化周波数： "+ Edit2->Text + " kHz");
    RichEdit1->Lines->Add("// 種類：　　　　 "+
                      RadioGroup1->Items->Strings[RadioGroup1->ItemIndex]);

    RichEdit1->Lines->Add("");
    buf = "//　　　　　　　　　　";
    for (int j=0; j<nBands; j++) buf = buf + "　　　　" + pLabelBnd[j]->Caption;
    RichEdit1->Lines->Add(buf);
    ParmOut("// 下側帯域端周波数 (kHz)", bands,   0, 2, fs);
    ParmOut("// 上側帯域端周波数 (kHz)", bands,   1, 2, fs);
    ParmOut("// ゲイン　　　　　　　　", gain,    0, 1);
    ParmOut("// 重み　　　　　　　　　", weights, 0, 1);
    ParmOut("// 偏差　　　　　　　　　", deviat,  0, 1);
    for (int j=0; j<nBands; j++) dev_dB[j] = 20.0*log10(gain[j] + deviat[j]);
    ParmOut("// 偏差 (dB)　　　　　　 ", dev_dB,  0, 1);

    CalcAmp(jtype, nTaps);      // 振幅特性の計算
    Button5Click(Sender);       // 指定された周波数における利得の表示
    PaintBox1->Invalidate();
    Button1->Caption = "設計";
}
//---------------------------------------------------------------------------

// “説明”ボタン，説明の画面を表示
void __fastcall TForm1::Button2Click(TObject */*Sender*/)
{
    Form2->ShowModal();
}
//---------------------------------------------------------------------------

// “結果の保存”ボタン
void __fastcall TForm1::Button3Click(TObject */*Sender*/)
{
    if (SaveDialog1->Execute() == true)
        RichEdit1->Lines->SaveToFile(SaveDialog1->FileName);
}
//---------------------------------------------------------------------------

// “終了”ボタン
void __fastcall TForm1::Button4Click(TObject */*Sender*/)
{
    Close();
}
//---------------------------------------------------------------------------

// “入力”ボタン（指定された周波数における利得の表示）
void __fastcall TForm1::Button5Click(TObject */*Sender*/)
{
    double fM, fr, sum, w;
    AnsiString buf;

    if (!calc_ok) return;               // まだ設計されていない
	if (!testDouble(Edit4, fr)) return;
    if (fr < 0.0 || fr > fs/2)
    {
        ShowMessage("周波数を 0 以上，標本化周波数／2 以下に設定してください．");
        Edit4->SetFocus();
        return;
    }

    fM = ((double)nTaps - 1)/2.0;
    w = 2.0*M_PI*fr/fs;
    sum = 0.0;
    if (jtype == 0) // 通過域／阻止域型
    {
        for (int n=0; n<=(nTaps/2-1); n++) sum = sum + hk[n]*cos(w*(fM-n));
        sum = (nTaps % 2) ? fabs(2.0*sum + hk[(nTaps - 1)/2]) : fabs(2.0*sum);
    }
    else            // 微分器，ヒルベルト変換器
    {
        for (int n=0; n<=(nTaps/2-1); n++) sum = sum + hk[n]*sin(w*(fM-n));
        sum = fabs(2.0*sum);
    }

    if (sum > 1.0E-30)
        Label10->Caption = buf.sprintf("%9.4f dB", 20.0*log10(sum));
    else
        Label10->Caption = "-∞ dB";
}
//---------------------------------------------------------------------------

// マウスの位置により Edit4 に周波数を入力し，dB 値を表示する
void __fastcall TForm1::PaintBox1MouseMove(TObject *Sender,
      TShiftState /*Shift*/, int X, int /*Y*/)
{
    AnsiString buf;

    if (!calc_ok) return;
    Edit4->SetFocus();
    if (X < X0d || X > X0d+5*Xfr) return;
    double val = fs*(X - X0d)/(10.0*Xfr);
    Edit4->Text = buf.sprintf("%-8.4f", val);
    Button5Click(Sender);
}
//---------------------------------------------------------------------------

// リターンキーで“入力”ボタンのクリックと同じ動作をさせる
void __fastcall TForm1::Edit4KeyPress(TObject *Sender, char &Key)
{
    if (Key == VK_RETURN)
    {
        Key = 0;
        Button5Click(Sender);
    }
}
//---------------------------------------------------------------------------

// 振幅特性描画の指令
void __fastcall TForm1::PaintBox1Paint(TObject */*Sender*/)
{
    if (calc_ok) DispSpec(PaintBox1->Canvas);
}
//---------------------------------------------------------------------------

// 座標軸の文字の描画指令
void __fastcall TForm1::PaintBox2Paint(TObject */*Sender*/)
{
    if (calc_ok) DispChr(PaintBox2->Canvas);
}
//---------------------------------------------------------------------------

// 標本化周波数変更に対するエディットボックス表示の処理
void __fastcall TForm1::Edit2Change(TObject */*Sender*/)
{
    if (Edit2->Text == "") return;
    pEditBdH[StrToInt(Edit3->Text)-1]->Text = FloatToStr(StrToFloat(Edit2->Text)/2);
    calc_ok = false;
}
//---------------------------------------------------------------------------

// 帯域の数変更に対するエディットボックス表示の処理
void __fastcall TForm1::Edit3Change(TObject */*Sender*/)
{
    bool enabled;

	if (!testInt(Edit3, nBands)) return;
    if ( (jtype == 0) && (nBands < 2) || (nBands > 5) )
    {
        ShowMessage("帯域の数は 2 以上，5 以下にしてください．");
        Edit3->SetFocus();
        return;
    }

    for (int j=0; j<MAX_BANDS; j++)
    {
        if (j<nBands) enabled = true;
        else          enabled = false;
        pEditBdL[j]->Enabled  = enabled;
        pEditBdH[j]->Enabled  = enabled;
        pEditGain[j]->Enabled = enabled;
        pEditWght[j]->Enabled = enabled;
        pLabelBnd[j]->Enabled = enabled;
    }
    calc_ok = false;
}
//---------------------------------------------------------------------------

// カーソルによる周波数，dB値の表示の On/Off
void __fastcall TForm1::LabelOnOff(TObject */*Sender*/)
{
    bool tf = GetLongHint(Application->Hint) == "PaintBox1" ? true : false;
    Label9->Visible  = tf;
    Label10->Visible = tf;
}
//---------------------------------------------------------------------------

// TEdit 配列の初期化
void __fastcall TForm1::SetEdit(TEdit *Ed[], int nBand, int x0, int y0, AnsiString s)
{
    for( int i=nBand-1; i>=0 ; i-- )
    {
        Ed[i] = new TEdit(this);
        Ed[i]->Parent = this;
        Ed[i]->Width = 68;
        Ed[i]->Height = 25;
        Ed[i]->Top =  y0;
        Ed[i]->Left = i*88 + x0;
        Ed[i]->Text = s;
        Ed[i]->TabOrder = 6;
        Ed[i]->OnChange = EditNChange;  // OnChange をEditNChange に割り当てる
        if (i<2) Ed[i]->Enabled = true;
        else     Ed[i]->Enabled = false;
    }
}
//---------------------------------------------------------------------------

// フィルタのタイプにより，帯域の数の変更を可能／不可能の設定
void __fastcall TForm1::RadioGroup1Click(TObject *Sender)
{
    bool tf;

    jtype = RadioGroup1->ItemIndex;
    if (jtype == 0) tf = true;
    else            tf = false;
    Label2->Enabled  = tf;
    Edit3->Enabled   = tf;
    UpDown1->Enabled = tf;
    Edit3->Text = tf ? 2 : 1;
    Edit3Change(Sender);
    if (jtype != 0)
    {
        pEditBdH[0]->Text = Edit2->Text/2;
        if (jtype == 2)     // Hilbert 変換器の場合
        {
            pEditGain[0]->Text = "1";
			pEditGain[0]->Enabled = false;
        }
        pEditWght[0]->Enabled = false;
    }
}
//---------------------------------------------------------------------------

// 振幅特性の縦軸の伸縮
void __fastcall TForm1::ComboBox1Change(TObject */*Sender*/)
{
    PaintBox1->Invalidate();
}
//---------------------------------------------------------------------------

// Edit の内容の変化に対する処理
void __fastcall TForm1::EditNChange(TObject */*Sender*/)    // ユーザ宣言
{
    calc_ok = false;
}
//---------------------------------------------------------------------------

// 設計パラメータの入力
bool __fastcall TForm1::InputSpec(VectorD &weights, VectorD &gain)
{
	if (!testInt(Edit1, nOrder)) return false;  // フィルタの次数
    if ( (nOrder > MAX_ORDER) || (nOrder < 2) )
    {
        ShowMessage("フィルタの次数を2 以上，" + IntToStr(MAX_ORDER) +
                    " 以下にしてください．");
        Edit1->SetFocus();
        return false;
    }
    nTaps = nOrder + 1;
	if (!testDouble(Edit2, fs)) return false;   // 標本化周波数
    jtype = RadioGroup1->ItemIndex;             // フィルタのタイプ設定
	if (!testInt(Edit3, nBands)) return false;  // 帯域数の読み込み
    for (int i=0; i<nBands; i++)
    {
		if (!testDouble(pEditBdL[i], bands[i*2])) return false;
		if (!testDouble(pEditBdH[i], bands[i*2+1])) return false;
		if (!testDouble(pEditGain[i], gain[i])) return false;
		if (!testDouble(pEditWght[i], weights[i])) return false;
    }
    for (int i=0; i<2*nBands; i++)
        if (bands[i] < 0.0 || bands[i] > fs/2.0)
        {
            ShowMessage("帯域端周波数を0以上，標本化周波数/2 以下にしてください．");
            if (i % 2) pEditBdH[i/2]->SetFocus();
            else       pEditBdL[i/2]->SetFocus();
            return false;
        }
    for (int i=1; i<2*nBands; i++)
        if (bands[i-1] > bands[i])
        {
            ShowMessage("帯域端周波数の設定が間違っています．");
            if (i % 2) pEditBdH[i/2]->SetFocus();
            else       pEditBdL[i/2]->SetFocus();
            return false;
        }
    if ((jtype == 1) && (gain[0] <= 0.0))       // 微分器の場合
    {
        ShowMessage("利得は 0 より大きな値を指定してください．");
        pEditGain[0]->SetFocus();
        return false;
    }
// 通過域／阻止域型で一番高い帯域の利得が 0 ではない奇数次のフィルタの設計を禁止
    if ((jtype == 0) && (gain[nBands-1] != 0.0) && (nOrder % 2))
    {
        ShowMessage("次数は偶数にしてください．");
        Edit1->SetFocus();
        return false;
    }
    Edit4->Text = pEditBdH[0]->Text;     // dB 値を表示する周波数を設定

    return true;
}
//---------------------------------------------------------------------------

// 振幅特性の計算
void __fastcall TForm1::CalcAmp(int jtype, int nTaps)
{
    double fM, sum, dw;
    int k, n, M;

    calc_ok = true;
    fM = ((double)nTaps - 1)/2.0;
    M = (nTaps - 1)/2;
    dw = M_PI/(double)fDiv;

    if (jtype == 0) // 通過域／阻止域型
        for (k=0; k<=fDiv; k++)
        {
            sum = 0.0;
            for (n=0; n<=(nTaps/2-1); n++) sum = sum + hk[n]*cos(k*dw*(fM-n));
            sum = (nTaps % 2) ? fabs(2.0*sum + hk[M]) : fabs(2.0*sum);
            LogA_k[k] = (sum < Amin) ? AmindB : 20.0*log10(sum);
        }
    else            // 微分器，ヒルベルト変換器
        for (k=0; k<=fDiv; k++)
        {
            sum = 0.0;
            for (n=0; n<=(nTaps/2-1); n++) sum = sum + hk[n]*sin(k*dw*(fM-n));
            sum = fabs(2.0*sum);
            LogA_k[k] = (sum < Amin) ? AmindB : 20.0*log10(sum);
        }
}
//---------------------------------------------------------------------------

// 振幅特性の表示
void __fastcall TForm1::DispSpec(TCanvas *Cv)
{
    double dx, dy;

    if (!calc_ok) return;

    dx = Xfr*5.0/(double)fDiv;
    dy = YdB/(abs(StrToInt(ComboBox1->Text))/8.0);

	Cv->Pen->Color = clGray;
// x軸
	for (int i=0; i<=5; i++) Line(Cv, X0d+i*Xfr, Y0d, X0d+i*Xfr, Y0d+8*YdB+1);
// y軸
	for (int i=0; i<=8; i++) Line(Cv, X0d, Y0d+i*YdB, X0d+5*Xfr, Y0d+i*YdB);
	DispChr(PaintBox2->Canvas);
// 振幅特性の描画
	Cv->Pen->Color = clBlue;
	Cv->Pen->Width = 2;
	Cv->MoveTo(X0d, Y0d-round(dy*LogA_k[0]));
	for (int i=1; i<=fDiv; i++)
		Cv->LineTo(X0d+round(i*dx), Y0d-round(dy*LogA_k[i]));

// 帯域端周波数の表示
	Cv->Pen->Color = clFuchsia;
	Cv->Pen->Width = 1;
    for (int j=0; j<nBands*2; j++)
        if (bands[j] > 0.0 && bands[j] < 0.5)
        {
            dx = round(bands[j]*10.0*Xfr);
            Line(Cv, X0d+dx, Y0d, X0d+dx, Y0d+8*YdB);
        }
}
//---------------------------------------------------------------------------

// 座標軸の文字の表示
void __fastcall TForm1::DispChr(TCanvas *Cv)
{
    AnsiString buf;

    double dBw = abs(StrToInt(ComboBox1->Text));
    Cv->Font->Name = "ＭＳ 明朝";
    for (int i=0; i<=5; i++)
        Cv->TextOut(X0d+i*Xfr-18, PaintBox1->Height+8,
                     buf.sprintf("%8.6f", (fs*i)/10.0));
    Cv->TextOut(X0d+230, PaintBox1->Height+24, "周波数 (kHz)");
    for (int i=0; i<=4; i++)
        Cv->TextOut(28, Y0d+2*i*YdB-6, buf.sprintf("%3.0f", -i*dBw/4.0));
    Cv->Font->Handle = CreateFont(-12, 0, 900, 0, 0, false, false, false,
                                      SHIFTJIS_CHARSET, 0, 0, 0, 0, NULL);
    Cv->TextOut(12, Y0d+YdB*4.5, "利得 (dB)");
    Cv->Font->Handle = CreateFont(-12, 0, 0, 0, 0, false, false, false,
                                      SHIFTJIS_CHARSET, 0, 0, 0, 0, NULL);
}
//---------------------------------------------------------------------------

// 設計パラメータ及びその他の結果の出力（Vector<double> の内容の出力）
void __fastcall TForm1::ParmOut(AnsiString as, const VectorD &x, int n1, int dn, double c)
{
    AnsiString buf = as;

    for (int j=0; j<nBands; j++) buf.cat_sprintf("%14.8f", c*x[j*dn+n1]);
    RichEdit1->Lines->Add(buf);
}
//---------------------------------------------------------------------------

// 係数表示部分の高さ，幅の変更
void __fastcall TForm1::FormResize(TObject */*Sender*/)
{
    RichEdit1->Width = width_of_RichEdit1 + (Form1->Width - width_of_Form1);
// 次のステートメントを入れない場合，高さを縮めた場合に縦のスクロールバーが現れない
    RichEdit1->Height = height_of_RichEdit1 + (Form1->Height - height_of_Form1) - 1;
    RichEdit1->Height = height_of_RichEdit1 + (Form1->Height - height_of_Form1);
}
//---------------------------------------------------------------------------

// ウィンドウの最小サイズの設定
void __fastcall TForm1::WMGetMinMaxInfo(TWMGetMinMaxInfo &Msg)
{
    Msg.MinMaxInfo->ptMinTrackSize.x = width_of_Form1;
    Msg.MinMaxInfo->ptMinTrackSize.y = height_of_Form1;
}

