//---------------------------------------------------------------------------
// 双一次ｚ変換による IIR フィルタの設計
//   設計可能なフィルタ：Butterworth, Chebyshev, Inverse Chebyshev, Elliptic
//
//      T. W. Parks and C. S. Burrus : “Digital Filter Design,”
//      Chapter 7, John Wiley & Sons, 1987 参照
//
// 設計結果の表示スペースをマウスにより拡大可能にしたバージョン
// 逆チェビシェフ特性の設計を追加
//
// Embarcadero C++ Builder 2010 で再ビルド	2010/02/23
//---------------------------------------------------------------------------
#include <complex.h>
#include <vcl.h>
#pragma hdrstop

#include "Unit1.h"
#include "Unit2.h"

#include "Elliptics.cpp"
#include "MK_misc.hpp"

using namespace std;

const int MAX_ORDER = 20;
const int MAX_COEF  = MAX_ORDER + 1;
const int A_SIZE    = MAX_ORDER*2 + 2;
const int fDiv      = 2048;     // 振幅特性描画の際の点数

complex<double> s_pz[A_SIZE];   // ｓ平面の極・零点
complex<double> z_pz[A_SIZE];   // ｚ平面の極・零点
double ak[MAX_COEF];            // 直接形に対応する伝達関数の分母の係数
double bk[MAX_COEF];            // 直接形に対応する伝達関数の分子の係数
double akm[2][MAX_COEF/2];      // 縦続形に対応する伝達関数の分母の係数
double bkm[3][MAX_COEF/2];      // 縦続形に対応する伝達関数の分子の係数
double gain;                    // 縦続形の利得定数
double scalec[MAX_COEF/2+1];    // スケーリング係数
double f_sample;                // 標本化周波数
double fcLow;                   // 遮断周波数（低域側）
double fcHigh;                  // 遮断周波数（高域側）
double dBrpl;                   // 通過域のリップル
double dBd;                     // 阻止域の減衰量
int    z_ord;                   // 伝達関数の次数
double LogA_k[fDiv+1];          // 計算された対数振幅特性の値

// 直線描画のインライン関数
void inline Line(TCanvas *Cv, int x1, int y1, int x2, int y2)
{ Cv->MoveTo(x1, y1); Cv->LineTo(x2, y2); }

double Modulus(double q);
void LambdaPlane(double order, double eps, double wr,
                 double &Kk, double &v0, double &k2, double &wc);
int sButChe(int order, int np, double a, double b);
int sButterworth(int order, int np);
int sChebyshev(int order, double phi, int np);
int sElliptic(int order, double m, double Kk, double v0, double k, double wc, int np);
void sPlane(int order, double phi, double k2, double Kk, double v0, double k, double wc,
            int &np, int &nz);
int Bilinear(int np, int nz, double cgam, double c, double cbp, double wc);
void AddZero(int jt);
void ToDirect(double cgam, double scale);
void Expand(double ck[], int order);
void ToCascade();
double AbsHw(double w);
void CalcAmp();
void Scaling();
double PeakPicking(int sct);
double Hw_n(double w, int sect);

//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
    : TForm(Owner)
{
    X0d =  50;                  // x軸の原点（周波数 = 0）
    Y0d =  24;                  // y軸の原点（0 dB）
    Xfr = 105;                  // x軸の1目盛り（(fs/2)/5）
    YdB =  58;                  // y軸の1目盛り（表示dBの幅/8）
    ComboBox1->ItemIndex = 4;   // 表示範囲を -80 dB までに設定
    EnableCtrlShape();
    EnableCtrlBand();

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

    Button1->Top = Form1->ClientHeight - 30;
    Button2->Top = Form1->ClientHeight - 30;
    Button3->Top = Form1->ClientHeight - 30;
    Button4->Top = Form1->ClientHeight - 30;
    Button6->Top = Form1->ClientHeight - 30;
    RichEdit1->Top = Form1->ClientHeight - 333;
    Panel1->Height = (Form1->ClientHeight - 39) < 541 ? 541 : Form1->ClientHeight - 39;

    calc_ok = false;
    scl1_ok = false;
    scl2_ok = false;
}
//---------------------------------------------------------------------------

// “設計実行”ボタン
void __fastcall TForm1::Button1Click(TObject *Sender)
{
    int    np, nz, order;
    double fa, delta, scale, ang, v0, k, k2, wc, wr, Kk, cbp, eps, cgam;

    if (!ParmInput(order, eps, scale, delta)) return;
    if (jBand == 1)     // HPF の場合
    {
        ang = fcHigh*M_PI/f_sample;
        fa = f_sample/2.0;
    }
    else
    {
        ang = (fcHigh - fcLow)*M_PI/f_sample;
        fa = fcHigh;
    }

    if (Shape != 3)             // Elliptic 以外の場合
    {
        // 周波数変換（LPF → BPF）
        wc = tan(ang);
        cgam = cos(M_PI*(fa + fcLow)/f_sample)/cos(ang);
        fa = 2.0*M_PI*fcHigh/f_sample;
        cbp = (cgam - cos(fa))/sin(fa);
    }
    else                        // Elliptic
        EllipInp(order, fa, eps, ang, cgam, Kk, k, v0, k2, wc, wr, cbp);

    sPlane(order, delta, k2, Kk, v0, k, wc, np, nz);
    AddZero(Bilinear(np, nz, cgam, tan(ang), cbp, wc));
    ToDirect(cgam, scale);
    XferFuncOut(RichEdit1->Lines);
    CalcAmp();
    Button5Click(Sender);       // 指定された周波数における利得の表示
    PaintBox1->Invalidate();
    scl1_ok = true;
    scl2_ok = true;
}
//---------------------------------------------------------------------------

// “スケーリング”ボタン
void __fastcall TForm1::Button2Click(TObject */*Sender*/)
{
    if (!scl1_ok)
    {
        ShowMessage("フィルタがまだ設計されていません．");
        return;
    }
    if (!scl2_ok)
    {
        ShowMessage("既にスケーリングは実行しました．");
        return;
    }
    Button2->Caption = "計算中";
    Button2->Update();
    Scaling();
    CascadeOut(RichEdit1->Lines, "スケーリングされた縦続形の係数", scalec[0]);
    Button2->Caption = "スケーリング";
    scl2_ok = false;
}
//---------------------------------------------------------------------------

// “説明”ボタン，説明の画面を表示
void __fastcall TForm1::Button3Click(TObject */*Sender*/)
{
    if (!FileExists("explain.rtf"))
    {
        ShowMessage("説明のファイル（explain.rtf）が見つかりません．");
        return;
    }
    Form2->ShowModal();
}
//---------------------------------------------------------------------------

// “結果の保存”ボタン
void __fastcall TForm1::Button4Click(TObject */*Sender*/)
{
    if (!calc_ok)
    {
        ShowMessage("フィルタがまだ設計されていません．");
        return;
    }
    if (SaveDialog1->Execute() == true)
        RichEdit1->Lines->SaveToFile(SaveDialog1->FileName);
}
//---------------------------------------------------------------------------

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

    if (!calc_ok) return;               // まだ設計されていない
	if (!testDouble(Edit7, fr)) return;
	fr = fr*1.0E3;
    if (fr < 0.0 || fr > f_sample/2)
    {
        ShowMessage("周波数を 0 以上，標本化周波数／2 以下に設定してください．");
        Edit7->SetFocus();
        return;
    }
    abs_hz = AbsHw(fr*2.0*M_PI/f_sample);
    if (abs_hz > 1.0E-30) Label13->Caption = buf.sprintf("%9.4f dB", 10.0*log10(abs_hz));
    else                  Label13->Caption = "-∞ dB";
}
//---------------------------------------------------------------------------

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

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

    if (!calc_ok) return;
    Edit7->SetFocus();
    if (X < X0d || X > X0d+5*Xfr) return;
    double val = 1.0E-3*f_sample*(X - X0d)/(10.0*Xfr);
    val = 0.01*((int)(val*100.0 + 0.5));    // 小数点２桁で四捨五入
    Edit7->Text = buf.sprintf("%-6.2f", val);
    Button5Click(Sender);
}
//---------------------------------------------------------------------------

// リターンキーで“入力”ボタンのクリックと同じ動作をさせる
void __fastcall TForm1::Edit7KeyPress(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::ComboBox1Change(TObject */*Sender*/)
{
    PaintBox1->Invalidate();
}
//---------------------------------------------------------------------------

// フィルタの形状により，入力可能な箇所を有効にする
void __fastcall TForm1::RadioGroup1Click(TObject */*Sender*/)
{
    EnableCtrlShape();
    calc_ok = false;
}
//---------------------------------------------------------------------------

// フィルタの通過／阻止帯域により，入力可能な箇所を有効にする
void __fastcall TForm1::RadioGroup2Click(TObject */*Sender*/)
{
    EnableCtrlBand();
    calc_ok = false;
}
//---------------------------------------------------------------------------

// 設計パラメータ変更に対応する処理
void __fastcall TForm1::OnChangeParm(TObject */*Sender*/)
{
    calc_ok = false;
}
//---------------------------------------------------------------------------

// RagioGroup1 による Label, Edit の有効／無効の制御
void __fastcall TForm1::EnableCtrlShape()
{
    switch (RadioGroup1->ItemIndex)
    {
        case 0: LabelEditTF(false, false); break;
        case 1: LabelEditTF( true, false); break;
        case 2: LabelEditTF(false, true ); break;
        case 3: LabelEditTF( true, true ); break;
    }
}
//---------------------------------------------------------------------------

// “通過域のリップル”，“阻止域の減衰量”の入力部の制御
void __fastcall TForm1::LabelEditTF(bool pass, bool reject)
{
    Label4->Enabled = pass;     // “通過域のリップル”
    Label5->Enabled = pass;
    Edit3->Enabled  = pass;
    Label6->Enabled = reject;   // “阻止域の減衰量”
    Label7->Enabled = reject;
    Edit4->Enabled  = reject;
}
//---------------------------------------------------------------------------

// RagioGroup2 による Label, Edit の有効／無効の制御
void __fastcall TForm1::EnableCtrlBand()
{
    bool tf = RadioGroup2->ItemIndex <= 1 ? false : true;
    Label10->Enabled = tf;
    Label11->Enabled = tf;
    Edit6->Enabled   = tf;
}
//---------------------------------------------------------------------------

// 設計パラメータの入力
bool __fastcall TForm1::ParmInput(int &order, double &eps, double &scale, double &delta)
{
    int min_order, max_order;

////    min_order = RadioGroup2->ItemIndex <= 1 ? 2 : 1;
    min_order = 1;
    max_order = RadioGroup2->ItemIndex <= 1 ? MAX_ORDER : MAX_ORDER/2;
	if (!testInt(Edit1, order)) return false;    // フィルタの次数
	if ( (order > max_order) || (order < min_order) )
	{
		ShowMessage("フィルタの次数を " + IntToStr(min_order) + " 以上，"
					+ IntToStr(max_order) + " 以下にしてください．");
		Edit1->SetFocus();
		return false;
	}
	if (!testDouble(Edit2, f_sample)) return false;     // 標本化周波数
	f_sample = fabs(f_sample)*1.0E3;
	Shape = RadioGroup1->ItemIndex;         // フィルタの形状
	jBand = RadioGroup2->ItemIndex;         // フィルタの通過域/阻止域による分類
	if (!testDouble(Edit5, fcHigh)) return false;       // 遮断周波数
	fcHigh = fcHigh*1000.0;                             // kHz 単位に換算
    if ((fcHigh <= 0.0) || (fcHigh >= f_sample/2.0))
    {
        ShowMessage("遮断周波数１を標本化周波数の 1/2 未満の正の値に設定してください．");
        Edit5->SetFocus();
        return false;
    }

    if (jBand >= 2)              // BPF, BRF の場合
    {
		if (!testDouble(Edit6, fcLow)) return false;    // もう一方の遮断周波数
        fcLow = fcLow*1000.0;                           // kHz 単位に換算
        if ((fcLow <= 0.0) || (fcLow >= f_sample/2.0))
        {
            ShowMessage("遮断周波数２を標本化周波数の 1/2 未満の正の値に設定してください．");
            Edit6->SetFocus();
            return false;
        }
    }
    else fcLow = 0.0;           // LPF, HPF の場合

    if (fcHigh < fcLow) swap(fcLow, fcHigh);
    Edit7->Text = FloatToStr(fcHigh/1.0e3);     // dB 値を表示する周波数を fcHigh に設定

    // 通過域のリップル
    if (Shape & 1)              // Chebyshev, Elliptic の場合
    {
		if (!testDouble(Edit3, dBrpl)) return false;    // 通過域のリップル
        if (dBrpl <= 0.0)
        {
            ShowMessage("通過域のリップルを 0 より大きな値に設定してください．");
            Edit3->SetFocus();
            return false;
        }
        scale = (order & 1) == 0 ? pow(10.0, dBrpl/20.0) : 1.0;
        if (Shape == 1) delta = pow(10.0, dBrpl/20.0);  // Chebyshev
        else eps = sqrt(pow(10.0, dBrpl/10.0) - 1.0);   // Elliptic
    }
    else
    {
    // 阻止域のリップル
        if (Shape == 2)
        {
            if (!dBdownInput()) return false;
            delta = pow(10.0, dBd/20.0);
        }
        scale = 1.0;
    }
    return true;
}
//---------------------------------------------------------------------------

// 阻止域の減衰量の入力
bool __fastcall TForm1::dBdownInput()
{
    if (!testDouble(Edit4, dBd)) return false;     // 阻止域の減衰量
    if ( dBd <= 0.0)
    {
        ShowMessage("阻止域の減衰量は正の値を入力してください．");
        Edit4->SetFocus();
        return false;
    }
    return true;
}
//---------------------------------------------------------------------------

// Elliptic フィルタの設計パラメータ入力その他
bool __fastcall TForm1::EllipInp(double order, double fa, double eps, double arg,
                                 double &cgam, double &Kk, double &k, double &u,
                                 double &k2, double &wc, double &wr, double &cbp)
{
    double ang, sang, cang, f3, tmp, m1, q, wrtan, wrtan2;

    cgam = cos((fa + fcLow)*M_PI/f_sample)/cos(arg);
    if (!dBdownInput()) return false;
    m1 = eps/sqrt(pow(10.0, dBd/10.0) - 1.0);
    m1 = m1*m1;
    k = Modulus(ellpC(m1)/(order*ellpC(1.0 - m1)));
    wr = (jBand & 1) ? k : 1.0/k;

    if (jBand <2)               // LPF, HPF の場合
        f3 = atan(tan(arg)*wr)*f_sample/M_PI;
    else                        // BPF, BRF の場合
    {
        wrtan = wr*tan(arg);
        wrtan2 = wrtan*wrtan;
        tmp = (cgam + sqrt(wrtan2*((1.0 - cgam*cgam) + wrtan2)))/(1.0 + wrtan2);
        f3 = (M_PI_2 - asin(tmp))*f_sample/(2.0*M_PI);
    }

    ang = f3*M_PI/f_sample;
    cang = cos(ang);
    sang = sin(ang);

    if (jBand <2)               // LPF, HPF の場合
        wr = sang/(cang*tan(arg));
    else                        // BPF, BRF の場合
    {
        q = cang*cang - sang*sang;
        sang = 2.0*cang*sang;
        cang = q;
        wr = (cgam - cang)/(sang*tan(arg));
    }
    if (jBand & 1) wr = 1.0/wr; // HPF, BRF の場合
    if (wr < 0.0) wr = -wr;
    cbp = wr;
    LambdaPlane(order, eps, wr, Kk, u, k2, wc);
    return true;
}
//---------------------------------------------------------------------------

// λ平面における位置を捜す
void LambdaPlane(double order, double eps, double wr,
                 double &Kk, double &v0, double &k2, double &wc)
{
    double  k, k1d;

    wc = 1.0;
    k = wc/wr;
    k2 = k*k;
    Kk = ellpC(1.0-k2);
    k1d = Modulus(order*ellpC(k2)/Kk);                          // k1'
    v0 = Kk*ellpI(1.0/eps, k1d)/(order*ellpC(1.0-k1d*k1d));     // 式(7.82)
}
//---------------------------------------------------------------------------

// sButterworth, sChebyshev で使用する関数
int sButChe(int order, int np, double a, double b)
{
    double theta;
    int nz = 0;

    theta = (order & 1) ? 0.0 : M_PI/(2.0*order);
    for (int j=0; j<np; j++)    // 極
    {
        s_pz[j] = complex<double>(-a*cos(theta), b*sin(theta));
        if (Form1->Shape == 2) s_pz[j] = 1.0/s_pz[j]; // s → 1/s
        theta = theta + M_PI/order;
    }
    if (Form1->Shape == 2)
    {
        for (int j=0; j<order; j++)
            s_pz[j+np] = complex<double>(0.0, 1.0/cos((2*j+1)*M_PI/(2.0*order)));
        nz = order/2;
    }
    if (Form1->jBand & 1)      // HPF, BRF の場合
    {
        for (int j=0; j<np; j++) s_pz[j] = 1.0/s_pz[j]; // s → 1/s
        nz = np;
        if (Form1->jBand == 3) nz = nz + order/2;
        if (Form1->Shape < 2)
            for (int j=0; j<nz; j++) s_pz[j+np] = 0.0;
        else
            for (int j=0; j<nz; j++) s_pz[j+np] = 1.0/s_pz[j+np];
    }
    return (nz);
}
//---------------------------------------------------------------------------

// Butterworth フィルタ：s 平面の極と零点の計算（遮断角周波数 = 1）
int sButterworth(int order, int np)
{
    return (sButChe(order, np, 1.0, 1.0));
}
//---------------------------------------------------------------------------

// Chebyshev フィルタ：s 平面の極と零点の計算（遮断角周波数 = 1）
int sChebyshev(int order, double delta, int np)
{
    double alpha, x;

    x = sqrt(delta*delta - 1.0);                // 逆Chebyshev
    if (Form1->Shape == 1) x = 1.0/x;           // Chebyshev
    alpha = log(x + sqrt(x*x  + 1.0))/order;    // arcsinh(x)
    return (sButChe(order, np, sinh(alpha), cosh(alpha)));
}
//---------------------------------------------------------------------------

// Elliptic フィルタ：s 平面の極と零点の計算（遮断角周波数 = 1）
int sElliptic(int order, double m, double Kk, double v0, double k, double wc, int np)
{
    int    nz;
    double denom, r, sn1, cn1, dn1, sn, cn, dn;

    nz = order/2;
    sncndn(v0, m, sn1, cn1, dn1);       // 式(7.85) sn', cn', dn' を求める
    for (int i=0; i<nz; i++)            // 零点
    {
        sncndn((order-1-2*i)*Kk/order, 1.0-m, sn, cn, dn);
        s_pz[i+np] = complex<double>(0.0, wc/(k*sn));
    }
    for (int i=0; i<np; i++)            // 極
    {
        sncndn((order-1-2*i)*Kk/order, 1.0-m, sn, cn, dn);
        r = k*sn*sn1;
        denom = cn1*cn1 + r*r;
        s_pz[i] = complex<double>(-wc*cn*dn*sn1*cn1/denom, wc*sn*dn1/denom);
    }
    if (Form1->jBand & 1)
    {
        for (int j=0; j<(np+nz); j++) s_pz[j] = 1.0/s_pz[j];
        while (np > nz) s_pz[np+nz++] = 0.0;
    }
    return (nz);
}
//---------------------------------------------------------------------------

// s 平面の極と零点の計算（遮断角周波数 = 1）
void sPlane(int order, double delta, double m, double Kk, double v0, double k, double wc,
            int &np, int &nz)
{
    for (int j=0; j<A_SIZE; j++) s_pz[j] = 0.0;
    np = (order+1)/2;
    switch (Form1->Shape)
    {
        case 0: nz = sButterworth(order, np);
                break;
        case 1:
        case 2:
                nz = sChebyshev(order, delta, np);
                break;
        case 3: nz = sElliptic(order, m, Kk, v0, k, wc, np);
                break;
    }
}
//---------------------------------------------------------------------------

// K/K' から母数（modulus）を求める
double Modulus(double u)
{
    double a, b, c, d, q;
    int j;

    q = exp(-M_PI*u);
    a = 1.0;
    b = 1.0;
    c = 1.0;
    d = q;
    for (j=0; j<15; j++)
    {
        a = a + 2.0*c*d;
        c = c*d*d;
        b = b + c;
        d = d*q;
        if (c < 1.0E-16) break;
    }
    if (j>=15) ShowMessage("Modulus: 収束しません．");
    a = b/a;
    return (4.0*sqrt(q)*a*a);
}
//---------------------------------------------------------------------------

// 双一次ｚ変換
int Bilinear(int np, int nz, double cgam, double c, double cbp, double wc)
{
    complex<double> r, root, cx, ca, cb, b4ac;
    double xc;
    int    nc, j, jt;

    xc = Form1->Shape == 3 ? c : wc;
    for (int i=0; i<A_SIZE; i++) z_pz[i] = 0.0;
    nc = np;
    jt = -1;
    j = 0;

    for (int i_pz=0; i_pz<2; i_pz++)    // s 平面から z 平面への変換
    {
        do
        {
            r = s_pz[j++];

            if (Form1->jBand < 2)       // LPF, HPF の場合
            {
                z_pz[++jt] = (1.0 + xc*r)/(1.0 - xc*r);
                if (imag(r) != 0.0)
                {
                    jt++;
                    z_pz[jt] = conj(z_pz[jt-1]);
                }
            }
            else                        // BPF, BRF の場合
            {
                cx = r*(Form1->Shape == 1 ? cbp : c);
                ca = 1.0 - cx;
                cb = -2.0*cgam;
                b4ac = sqrt(cb*cb - 4.0*(1.0 - cx*cx));
                root = (b4ac - cb)/(2.0*ca);

                z_pz[++jt] = root;
                if (imag(root) != 0.0)
                    z_pz[++jt] = conj(root);
                if ((imag(r) != 0.0) || (imag(root) == 0))
                {
                    root = -(b4ac + cb)/(2.0*ca);
                    z_pz[++jt] = root;
                    if (imag(root) != 0.0) z_pz[++jt] = conj(root);
                }
            }
        } while (--nc > 0);

        if (i_pz == 0)
        {
            z_ord = jt+1;
            if (nz <= 0)
            {
                if (Form1->Shape <= 1) return (jt);
                else                   break;
            }
        }
        nc = nz;
    }
    return (jt);
}
//---------------------------------------------------------------------------

// 残りの零点を付加
void AddZero(int jt)
{
    if (Form1->Shape <= 2)          // Butterworth，Chebyshev
        while (2*z_ord - 1 > jt)
        {
            if (Form1->jBand != 1) z_pz[++jt] = -1.0;
            if (Form1->jBand == 1 || Form1->jBand == 2) z_pz[++jt] = 1.0;
        }
    else                            // Elliptic
        while (2*z_ord - 1 > jt)
        {
            z_pz[++jt] = -1.0;
            if (Form1->jBand >= 2) z_pz[++jt] = 1.0;
        }
}
//---------------------------------------------------------------------------

// 直接形の係数に変換
void ToDirect(double cgam, double scale)
{
    double a, ai, an, bn, gam;
    int    mh, jh, jl;

    Expand(ak, 0);
    Expand(bk, z_ord);
    a = 1.0;

    if (Form1->jBand == 2)      // BPF の場合
    {
        gam = M_PI_2 - asin(cgam);
        mh = z_ord/2;
        an = ak[mh];
        bn = bk[mh];
        ai = 0.0;
        if (mh > ((z_ord/4)*2))
        {
            ai = 1.0;
            an = 0.0;
            bn = 0.0;
        }
        for (int j=1; j<=mh; j++)
        {
            a = gam*j - ai*M_PI_2;
            jh = mh + j;
            jl = mh - j;
            an = an + cos(a)*(ak[jh] + (1.0 - 2.0*ai)*ak[jl]);
            bn = bn + cos(a)*(bk[jh] + (1.0 - 2.0*ai)*bk[jl]);
        }
    }
    else                        // BPF 以外
    {
        if (Form1->jBand == 1) a = -1.0;
        an = 1.0;
        bn = 1.0;
        for (int j=1; j<=z_ord; j++)
        {
            an = a*an + ak[j];
            bn = a*bn + bk[j];
        }
    }

    gain = an/(bn*scale);
    if ((Form1->Shape != 3) && (bn == 0)) gain = 1.0;

    for (int j=0; j<=z_ord; j++) bk[j] = gain*bk[j];
}
//---------------------------------------------------------------------------

// 多項式の係数に展開
void Expand(double ck[], int order)
{
    complex<double> xk[A_SIZE];

    xk[0] = 1.0;
    for (int j=1; j<MAX_COEF; j++) xk[j] = 0.0;
    for (int j=0; j<z_ord; j++)
        for (int i=0; i<=j; i++)
            xk[j-i+1] = xk[j-i+1] - z_pz[j+order]*xk[j-i];
    for (int j=0; j<=z_ord; j++) ck[j] = real(xk[j]);
}
//---------------------------------------------------------------------------

// 縦続形の係数に変換
void ToCascade()
{
    int m, jz, pos;
    double min, tmp_a0, tmp_a1, tmp_b0, tmp_b1, tmp_b2;

    m = 0;
    for (int j=0; j<z_ord; j++)
    {
        jz = j + z_ord;
        if (fabs(imag(z_pz[j])) > 1.0E-16)    // 複素極のステージ
        {
            akm[0][m] = real(2.0*z_pz[j]);              // a1m
            akm[1][m] = -norm(z_pz[j]);                 // a2m
            bkm[0][m] = 1.0;                            // b0m
            bkm[1][m] = -real(z_pz[jz] + z_pz[jz+1]);   // b1m
            bkm[2][m] = real(z_pz[jz]*z_pz[jz+1]);      // b2m
            j++;
        }
        else                                // 実極のステージ
        {
            if (z_ord%2 != 0)               // 奇数次の場合
            {
                akm[0][m] = real(z_pz[j]);                  // a1m
                akm[1][m] = 0.0;                            // a2m = 0
                bkm[0][m] = 1.0;                            // b0m
                bkm[1][m] = real(-z_pz[jz]);                // b1m
                bkm[2][m] = 0.0;                            // b2m = 0
            }

            else                            // 偶数次の場合
            {
                akm[0][m] = real(z_pz[j] + z_pz[j+1]);      // a1m
                akm[1][m] = -real(z_pz[j]*z_pz[j+1]);       // a2m
                bkm[0][m] = 1.0;                            // b0m
                bkm[1][m] = -real(z_pz[jz] + z_pz[jz+1]);   // b1m
                bkm[2][m] = real(z_pz[jz]*z_pz[jz+1]);      // b2m
                j++;
            }

        }
        m++;
    }
// ソート，極の絶対値が小さい順，ただし実極は先頭になる
    for (int i=0; i<(z_ord+1)/2-1; i++)
    {
        pos = i;
        min = fabs(akm[1][i]);
        for (int j=i; j<(z_ord+1)/2; j++)
        {
            if (fabs(akm[1][j]) <= min)
            {
                pos = j;
                min = fabs(akm[1][j]);
                tmp_a0 = akm[0][j];
                tmp_a1 = akm[1][j];
                tmp_b0 = bkm[0][j];
                tmp_b1 = bkm[1][j];
                tmp_b2 = bkm[2][j];
            }
        }
        akm[0][pos] = akm[0][i];
        akm[1][pos] = akm[1][i];
        bkm[0][pos] = bkm[0][i];
        bkm[1][pos] = bkm[1][i];
        bkm[2][pos] = bkm[2][i];
        akm[0][i] = tmp_a0;
        akm[1][i] = tmp_a1;
        bkm[0][i] = tmp_b0;
        bkm[1][i] = tmp_b1;
        bkm[2][i] = tmp_b2;
    }
}
//---------------------------------------------------------------------------

// 伝達関数の係数などの出力
void __fastcall TForm1::XferFuncOut(TStrings *ML)
{
    AnsiString buf;

    ML->Clear();
    switch (RadioGroup2->ItemIndex)
    {
        case 0: ML->Add("低域通過フィルタ"); break;
        case 1: ML->Add("高域通過フィルタ"); break;
        case 2: ML->Add("帯域通過フィルタ"); break;
        case 3: ML->Add("帯域阻止フィルタ"); break;
    }
    ML->Add(RadioGroup1->Items->Strings[RadioGroup1->ItemIndex] + "特性");
    ML->Add(buf.sprintf("次数　　　　：%d 次", z_ord));
    ML->Add(buf.sprintf("標本化周波数：%8.6f kHz", f_sample/1.0E3));
    if (RadioGroup2->ItemIndex < 2)
        ML->Add(buf.sprintf("遮断周波数　：%8.6f kHz", fcHigh/1.0E3));
    else
    {
        ML->Add(buf.sprintf("遮断周波数１：%8.6f kHz", fcLow/1.0E3));
        ML->Add(buf.sprintf("遮断周波数２：%8.6f kHz", fcHigh/1.0E3));
    }
    if ((RadioGroup1->ItemIndex & 1) == 1)
        ML->Add(buf.sprintf("通過域のリップル：%5.2f dB", dBrpl));
    if (RadioGroup1->ItemIndex >= 2)
        ML->Add(buf.sprintf("阻止域の減衰量　：%5.2f dB", dBd));

    ML->Add("");
    ML->Add("直接形の係数");
    ML->Add(" m    分母の係数 (am)    分子の係数 (bm)");
    // 直接形の係数の出力
    ML->Add(buf.sprintf("%2d  %36.10E", 0, bk[0]));
    for (int j=1; j<=z_ord; j++)
        ML->Add(buf.sprintf("%2d %18.10E %18.10E", j, -ak[j], bk[j]));

    // 縦続形の係数の出力
    ToCascade();
    CascadeOut(RichEdit1->Lines, "縦続形の係数", 0.0);
    ML->Add("");
    ML->Add(buf.sprintf("利得定数 = %16.10E", gain));
}
//---------------------------------------------------------------------------

// 縦続形の係数の出力
void __fastcall TForm1::CascadeOut(TStrings *ML, AnsiString s, double scl)
{
    AnsiString buf;

    ML->Add("");
    ML->Add(s);
    if (scl != 0.0)
    {
        ML->Add("");
        ML->Add(buf.sprintf("s0  = %17.10E", scl));
    }
    for (int j=0; j<(z_ord+1)/2; j++)
    {
        ML->Add("");
        ML->Add(buf.sprintf("a1%d = %17.10E", j+1, akm[0][j]));
        if (akm[1][j] != 0.0)       // 複素極のステージ
            ML->Add(buf.sprintf("a2%d = %17.10E", j+1, akm[1][j]));

        ML->Add(buf.sprintf("b0%d = %17.10E", j+1, bkm[0][j]));
        ML->Add(buf.sprintf("b1%d = %17.10E", j+1, bkm[1][j]));
        if (akm[1][j] != 0.0)       // 複素極のステージ
            ML->Add(buf.sprintf("b2%d = %17.10E", j+1, bkm[2][j]));
    }
}
//---------------------------------------------------------------------------

// 正規化角周波数に対する振幅特性の値
double AbsHw(double w)
{
    complex<double> z_inv, numer, denom, uj(0,-1);
    z_inv = exp(uj*w);
    denom = ak[z_ord];
    for (int n=z_ord-1; n>=0; n--) denom = denom*z_inv + ak[n];
    numer = bk[z_ord];
    for (int n=z_ord-1; n>=0; n--) numer = numer*z_inv + bk[n];
    return (norm(numer/denom));
}
//---------------------------------------------------------------------------

// 振幅特性の計算
void CalcAmp()
{
    double norm_hz, dw;

    Form1->calc_ok = true;
    dw = M_PI/fDiv;
    for (int k=0; k<=fDiv; k++)
    {
        norm_hz = AbsHw(k*dw);
        LogA_k[k] = norm_hz < 1.0E-10 ? -100 : 10.0*log10(norm_hz);
    }
}
//---------------------------------------------------------------------------

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

    if (!calc_ok) return;

    dx = Xfr*5.0/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;
    DispFc(PaintBox1->Canvas, fcLow, Xfr, X0d, Y0d, YdB);
    DispFc(PaintBox1->Canvas, fcHigh, Xfr, X0d, Y0d, YdB);
}
//---------------------------------------------------------------------------

// 遮断周波数の表示
void __fastcall TForm1::DispFc(TCanvas *Cv, double fc, int Xfr, int X0d, int Y0d, int YdB)
{
    if (fc <= 0.0 || fc >= f_sample/2.0) return;
    double x = round(10.0*Xfr*fc/f_sample + X0d);
    Line(Cv, x, Y0d, x, 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, Y0d+YdB*8+8,
                     buf.sprintf("%8.6f", (f_sample*i)/1.0E4));
    Cv->TextOut(X0d+230, Y0d+YdB*8+26, "周波数 (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);
}
//---------------------------------------------------------------------------

// L∞ノルムによるスケーリング
void Scaling()
{
    double norm[MAX_COEF/2];

    for (int j=0; j<(z_ord+1)/2; j++) norm[j] = PeakPicking(j);

    scalec[0] = 1.0/norm[0];
    for (int j=1; j<(z_ord+1)/2; j++) scalec[j] = norm[j-1]/norm[j];
    scalec[(z_ord+1)/2] = gain*norm[(z_ord+1)/2-1];

    for (int j =0; j<(z_ord+1)/2; j++)
        for (int k=0; k<3; k++)
            bkm[k][j] = scalec[j+1]*bkm[k][j];
}
//---------------------------------------------------------------------------

// 振幅特性の最大値検出
double PeakPicking(int sct)
{
    double y, w, amp_max, dw;

    amp_max = 0.0;
    dw = M_PI/1.E4;

    for (w=0; w<=M_PI; w+=dw)
    {
            y = Hw_n(w, sct);
        amp_max = amp_max < y ? y : amp_max;
    }
    return amp_max;
}
//---------------------------------------------------------------------------

// 入力点からオーバーフローに注目する点までの伝達関数
double Hw_n(double w, int sect)
{
    complex<double> h, z, z2, uj(0, -1);

    z =exp(uj*w);
    z2 = z*z;
    h = 1.0/(1.0 - akm[0][0]*z - akm[1][0]*z2); // 第１セクション
    if (sect==0) return ( abs(h) );
    for (int j=1; j<=sect; j++)                 // 第２セクション以降
        h = h*(bkm[0][j-1] + bkm[1][j-1]*z + bkm[2][j-1]*z2)
             /(1.0 - akm[0][j]*z - akm[1][j]*z2);
    return ( abs(h) );
}
//---------------------------------------------------------------------------


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;
}

