極・零点の数値計算のソースコード

  • 2021年6月28日
  • Home

ファイル構成

ファイルの中身

4_zero_pole_calculation.cpp

#include <fstream>
#include <string>
#include <vector>
#include <iostream>
#include <complex>
#include "csv_import.h"
#include "4_zero_pole_calculation_function.h"

using namespace std;

// =============================================================================================================================================
// 実行プログラム
// =============================================================================================================================================
int main()
{
    //------------------------------------------------------------------------------------------------------------------------------------------
    // インターフェース
    string arx_model_coefficient_file_name;
    cout << "Enter arx model coefficient data file name" << endl;
    cin  >> arx_model_coefficient_file_name;
    
    //------------------------------------------------------------------------------------------------------------------------------------------
    // arxモデルの係数取り込み
    ifstream inputfile( "./data_strage/" + arx_model_coefficient_file_name + ".csv" );
    string   line;
    vector<string> strvec;
    getline( inputfile, line ); // 第1行目を捨てる

    getline( inputfile, line ); 
    strvec = split(line, ',');
    int output_order = stof( strvec[0] );
    int  input_order = stof( strvec[1] );
    int    dead_time = stof( strvec[2] );

    vector<float> output_coef;
    output_coef.push_back( 1 );
    for( int i=0; i<output_order; i++ )
    {
        // 出力の係数
        output_coef.push_back( stof( strvec[i+4] ) );
    }

    vector<float> input_coef;
    for( int i=0; i<input_order; i++ )
    {
        // 入力の係数
        input_coef.push_back( stof( strvec[i+4+output_order] ) );
    }
    for( int i=0; i<dead_time; i++ )
    {
        input_coef.push_back( 0 );
    }

    //--------------------------------------------------------------------------------------------------------------------------------------------
    // 極の導出
    vector< complex<float> > pole( output_order );
    pole = polynomial_roots_calculation( output_order, output_coef );

    //--------------------------------------------------------------------------------------------------------------------------------------------
    // 零点の導出
    vector< complex<float> > zero( input_order );
    zero = polynomial_roots_calculation( input_order, input_coef );

    // --------------------------------------------------------------------------------------------------------------------------------------------
    // CSV出力(arxモデルの極とゼロを出力)
    ofstream outputfile("./data_strage/4_poles_and_zeros.csv");

    // 第1行目の記入
    outputfile << "poles list" << "," << "," << "," << "zeros list" << endl;

    // 第2行目の記入
    outputfile << "Real part" << "," << "Imaginary part" << "," << ",";
    outputfile << "Real part" << "," << "Imaginary part" << endl ;
    
    // 極と零点の出力
    for ( int i=0; i<10; i++ )
    {
        // 極の出力
        if( i<pole.size() )
        {
            outputfile << real( pole[i] ) << "," << imag( pole[i] ) << "," << ",";
        }
        else
        {
        } 

        // 零点の出力
        if( i<zero.size() )
        {
            outputfile << real( zero[i] ) << "," << imag( zero[i] );
        }
        else
        {
        }
        outputfile << endl;
        
    }
    cout << "poles and zeros calculation is completed" << endl;
}

4_zero_pole_calculation_function.h

#include <vector>
#include <complex>
using namespace std;

// 多項式の根の計算
vector< complex<float> > polynomial_roots_calculation( int order, vector<float> coefficient );

// 多項式の値
complex<float> polynomial_value( int order, vector<float> coefficient, complex<float> value );

// 多項式の値
complex<float> polynomial_derivative( int order, vector<float> coefficient, complex<float> value );

polynomial_roots_calculation.cpp

#include <vector>
#include <complex>
#include <bits/stdc++.h>

#include "4_zero_pole_calculation_function.h"

using namespace std;

const double  pi = acos(-1);

//---------------------------------------------------------------------------------------
// 多項式の根の計算
vector< complex<float> > polynomial_roots_calculation( int order, vector<float> coefficient )
{
    vector< complex<float> > result( order ); // 多項式の根
    complex<float> value;                     // 多項式の値
    complex<float> derivative;                // 多項式の導関数

    // 初期条件のパラメータ
    float a1 = coefficient[1];
    float a0 = coefficient[0];
    float  R = 100;
    int    n = order;

    // 多項式の根の導出
    for ( int root_number=0; root_number<order; root_number++ )
    {
        // 初期値の計算
        float initial_real = -a1/(n*a0) + R*cos( 2*pi/n*(root_number-3/4) );
        float initial_imag =              R*sin( 2*pi/n*(root_number-3/4) );
        complex<float> initial( initial_real, initial_imag );
        
        // 初期値に対する多項式の値
        result[root_number] = initial;
        value = polynomial_value( order, coefficient, result[root_number] );

        // 多項式の1つの根を求めるループ
        while( abs(value)>0.001 )
        {
            value      =      polynomial_value( order, coefficient, result[root_number] ); // result[root_number]に対する多項式の値
            derivative = polynomial_derivative( order, coefficient, result[root_number] ); // result[root_number]に対する導関数

            result[root_number] = result[root_number] - value/derivative; // ニュートン法の式
        }
    }

    return result;
}

polynomial_value.cpp

#include <vector>
#include <complex>
#include <bits/stdc++.h>
using namespace std;

//---------------------------------------------------------------------------------------
// ある値での多項式の値の計算
complex<float> polynomial_value( int order, vector<float> coefficient, complex<float> value )
{
    int term_size = coefficient.size(); // 多項式の項の数
    complex<float> result; // ある値での多項式の値
    complex<float> term; // 各項の値

    // ある値での多項式の値を求めるループ
    for ( int term_number=0; term_number<term_size; term_number++ )
    {
        term = term + coefficient[term_number] * pow( value, order-term_number ); 
    }
    result = term;

    return result;
}

polynomial_derivative.cpp

#include <vector>
#include <complex>
#include <bits/stdc++.h>
using namespace std;

//---------------------------------------------------------------------------------------
// ある値での多項式の導関数の導出
complex<float> polynomial_derivative( int order, vector<float> coefficient, complex<float> value )
{
    int term_size = coefficient.size() - 1; // 多項式の項の数
    complex<float> result; // ある値での多項式の導関数
    complex<float> term; // 各項の値

    // ある値での多項式の導関数を求めるループ
    for ( int term_number=0; term_number<term_size; term_number++ )
    {
        term = term + coefficient[term_number] * ( order-term_number ) * pow( value, order-term_number-1 ); 
    }
    result = term;

    return result;
}