#include <iostream>
#include <math.h>
using namespace std;
double numerador[8] = { 0 };
double denominador[8] = { 0 };
double Binomial_E (double a, int n)
{
if (n == 0.0)
return 1.0;
else
return (1.0 - (1.0 + a) / n) * Binomial_E (a, n - 1.0);
}
double Binomial (double a, int n)
{
if (n == 0.0)
return 1.0;
else
return ((a + 1.0) / n - 1.0) * Binomial (a, n - 1.0);
}
void IIR_Pade (double alpha, double T, int orderPade, double operador, double A[],
double B[])
{
int L = (2 * orderPade) + 1;
double h[15] = { 0 };
double Gain = pow ((1 / T), alpha);
for (int n = 0; n <= L; n++)
{
for (int k = 0; k <= n; k++)
{
h[n] +=pow (operador, -alpha - n + k) * pow (1.0 - operador,n - k) * Binomial_E (alpha,k) *Binomial (-alpha, n - k);
}
h[n]=h[n]*Gain;
}
int LM = 2 * orderPade;
double H2[orderPade][LM];
double temp;
for (int i = 0; i < orderPade; i++)
{
for (int j = 0; j < LM; j++)
{
if (j == (i + orderPade))
{
H2[i][j] = 1;
}
else
{
H2[i][j] = 0;
}
if (j < orderPade)
{
H2[i][j] = h[i - j + orderPade];
}
}
}
for (int i = orderPade - 1; i > 0; i--)
{
if (H2[i - 1][0] < H2[i][0])
for (int j = 0; j < LM; j++)
{
temp = H2[i][j];
H2[i][j] = H2[i - 1][j];
H2[i - 1][j] = temp;
}
}
for (int i = 0; i < orderPade; i++)
{
for (int j = 0; j < orderPade; j++)
{
if (j != i)
{
temp = H2[j][i] / H2[i][i];
for (int k = 0; k < LM; k++)
{
H2[j][k] -= H2[i][k] * temp;
}
}
}
}
for (int i = 0; i < orderPade; i++)
{
temp = H2[i][i];
for (int j = 0; j < LM; j++)
{
H2[i][j] = H2[i][j] / temp;
}
}
B[0] = 1;
for (int i = 0; i < orderPade; i++)
{
for (int j = 0; j < orderPade; j++)
{
B[i + 1] -= H2[i][j + orderPade] * h[j + orderPade + 1];
}
}
double H1[orderPade + 1][orderPade + 1];
for (int i = 0; i <= orderPade; i++)
{
for (int j = 0; j <= orderPade; j++)
{
if (j <= i)
{
H1[i][j] = h[i - j];
}
else
{
H1[i][j] = 0;
}
}
}
for (int i = 0; i <= orderPade; i++)
{
for (int j = 0; j <= orderPade; j++)
{
A[i] += H1[i][j] * B[j];
}
cout << A[i] << " ";
}
cout << endl;
for (int i = 0; i <= orderPade; i++)
{
cout << B[i] << " ";
}
}
int main ()
{
IIR_Pade (1, 0.001, 1, 1, numerador, denominador);
return 0;
}
