code. compile. run. debug. share.
#include "fractionalOrderFilter.h"
#include <fstream>
#include <iomanip>
// Con formato de matricial de Mathematica
void _guardarVector(const std::string& filename, const double* data, int size) {
std::ofstream file(filename);
file << "{";
for (int i = 0; i < size - 1; ++i) {
file << std::setprecision(15) << data[i] << "," << std::endl;
}
file << std::setprecision(15) << data[size - 1] << "}";
}
// Función para guardar los coeficientes de numerador y denominador en un archivo
// facilmente usados en Mathematica
void _guardarCoeficientes(const std::string& filename, const double* num, int ordenNum, const double* den, int ordenDen) {
std::ofstream file(filename);
file << "{{";
for (int i = 0; i < ordenNum; ++i) {
file << std::setprecision(15) << num[i] << "," << std::endl;
}
file << std::setprecision(15) << num[ordenNum] << "},{";
for (int i = 0; i < ordenDen; ++i) {
file << std::setprecision(15) << den[i] << "," <<std::endl;
}
file << std::setprecision(15) << den[ordenDen] << "}}";
}
int main() {
// Inicialización del filtro con parámetros comunes
FractionalOrderFilter filter(-0.5, 0.01, 0.5);
// Ejemplo de coeficientes de la respuesta impulso
int cantidad1 = 100;
double h[cantidad1];
filter._respuestaImpulso(cantidad1, h);
_guardarVector("datos1.txt", h, cantidad1);
// Ejemplo de coeficientes de Padé
int ordenNum1 = 9, ordenDen1 = 9;
double num1[ordenNum1 + 1], den1[ordenDen1 + 1];
filter._Pade(num1, ordenNum1, den1, ordenDen1);
_guardarCoeficientes("datos2.txt", num1, ordenNum1, den1, ordenDen1);
// Ejemplo de coeficientes de Shank
int ordenNum3 = 9, ordenDen3 = 9, cantidad3 = 1000;
double num3[ordenNum3 + 1], den3[ordenDen3 + 1];
filter._Shank(num3, ordenNum3, den3, ordenDen3, cantidad3);
_guardarCoeficientes("datos3.txt", num3, ordenNum3, den3, ordenDen3);
// Ejemplo de coeficientes de Prony
int ordenNum2 = 9, ordenDen2 = 9, cantidad2 = 1000;
double num2[ordenNum2 + 1], den2[ordenDen2 + 1];
filter._Prony(num2, ordenNum2, den2, ordenDen2, cantidad2);
_guardarCoeficientes("datos4.txt", num2, ordenNum2, den2, ordenDen2);
return 0;
}
#ifndef FRACTIONAL_ORDER_FILTER_H
#define FRACTIONAL_ORDER_FILTER_H
#include <vector>
#include <string>
const double _MIN_ORDER = -1.0;
const double _MAX_ORDER = 1.0;
const double _MIN_TYPE = 0.0;
const double _MAX_TYPE = 1.0;
struct _GlobalConfig {
double _orden;
double _periodo;
double _tipo;
};
struct _Scale {
double _valor;
bool _aplicar;
};
struct _MatrixParams {
double* _matriz;
double* _h;
double* _g_o_h;
int _filas;
int _inicio;
int _fin;
double* _result;
bool _positivo;
};
class FractionalOrderFilter {
public:
FractionalOrderFilter(double orden, double periodo, double tipo);
void _respuestaImpulso(int cantidad, double* respuesta);
void _Pade(double* num, int ordenNum, double* den, int ordenDen);
void _Shank(double* num, int ordenNum, double* den, int ordenDen, int cantidad);
void _Prony(double* num, int ordenNum, double* den, int ordenDen, int cantidad);
private:
_GlobalConfig _gC;
_Scale _Escala;
void _numPade(double* h, double* num,double ordenNum,double* den,double ordenDen);
void _Pade(double* h, double* num, int ordenNum, double* den, int ordenDen);
void _Shank(double* h, double* num, int ordenNum, double* den, int ordenDen, int cantidad);
void _Prony(double* h, double* num, int ordenNum, double* den, int ordenDen, int cantidad);
void _gaussJordan(double* A, double* b, int N);
void _swapRows(double* A, double* b, int N, int row1, int row2);
void _calcularMatriz(const _MatrixParams& params);
void _ponerCeros(double* array, int size);
};
#endif
#include "fractionalOrderFilter.h"
#include <cmath>
#include <algorithm>
FractionalOrderFilter::FractionalOrderFilter(double orden, double periodo, double tipo) {
if (orden < _MIN_ORDER || orden > _MAX_ORDER || tipo < _MIN_TYPE || tipo > _MAX_TYPE) {
return;
}
_gC._orden = orden;
_gC._periodo = periodo;
_gC._tipo = tipo;
_Escala._valor = pow(_gC._periodo, -_gC._orden);
_Escala._aplicar = true;
}
void FractionalOrderFilter::_swapRows(double* A, double* b, int N, int row1, int row2) {
if (row1 == row2) return;
for (int i = 0; i < N; i++) {
double temp = A[row1 * N + i];
A[row1 * N + i] = A[row2 * N + i];
A[row2 * N + i] = temp;
}
double temp = b[row1];
b[row1] = b[row2];
b[row2] = temp;
}
void FractionalOrderFilter::_gaussJordan(double* A, double* b, int N) {
for (int i = 0; i < N; i++) {
// Pivoteo parcial
int maxRow = i;
for (int k = i + 1; k < N; k++) {
if (fabs(A[k * N + i]) > fabs(A[maxRow * N + i])) {
maxRow = k;
}
}
_swapRows(A, b, N, i, maxRow);
double pivot = A[i * N + i];
for (int j = i; j < N; j++) {
A[i * N + j] /= pivot;
}
b[i] /= pivot;
for (int j = 0; j < N; j++) {
if (j != i) {
double factor = A[j * N + i];
for (int k = i; k < N; k++) {
A[j * N + k] -= factor * A[i * N + k];
}
b[j] -= factor * b[i];
}
}
}
}
void FractionalOrderFilter::_calcularMatriz(const _MatrixParams& params) {
for (int i = 0; i < params._filas; i++) {
for (int j = 0; j <= i; j++) {
double sumMatrix = 0.0;
double sumResult = 0.0;
for (int k = params._inicio; k < params._fin; k++) {
double val_i = (k - i >= 0) ? params._g_o_h[k - i] : 0.0;
double val_j = (k - j >= 0) ? params._g_o_h[k - j] : 0.0;
sumMatrix += val_i * val_j;
if (j == 0) {
sumResult += val_i * (params._positivo ? params._h[k] : params._h[k + 1]);
}
}
params._matriz[params._filas * i + j] = sumMatrix;
if (i != j) {
params._matriz[params._filas * j + i] = sumMatrix;
}
if (j == 0) {
params._result[i] = (params._positivo ? sumResult : -sumResult);
}
}
}
}
void FractionalOrderFilter::_ponerCeros(double* array, int size){
std::fill(array, array + size, 0.0);
}
void FractionalOrderFilter::_respuestaImpulso(int cantidad, double* respuesta) {
if (cantidad <= 0) {
return;
}
double k1 = pow(_gC._tipo, -_gC._orden);
double k2 = (1.0 - _gC._tipo) / _gC._tipo;
double inv_k2 = 1/k2;
double k2_power = 1.0;
double c1 = 1 + _gC._orden;
double c2 = 1 - _gC._orden;
double EE1 = 1.0;
double DD[cantidad];
double EE[cantidad];
DD[0] = 1.0;
EE[0] = 1.0;
for (int n = 0; n < cantidad; n++) {
if (n > 0){
EE[n] = (1 - c1 / n) * EE[n - 1];
EE1 *= (1 - c2 / n);
DD[n] = EE1 * (1 - 2 * (n & 1));
k2_power *= k2;
}
double sum = 0.0;
if (_gC._tipo == 1){
sum = EE[n];
}else{
double k2_local = k2_power;
for (int k = 0; k <= n; k++) {
sum += k2_local * EE[k] * DD[n - k];
k2_local *=inv_k2;
}
}
respuesta[n] =k1*sum*((_Escala._aplicar == true) ? _Escala._valor : 1);
}
}
void FractionalOrderFilter::_Pade(double* num, int ordenNum, double* den, int ordenDen) {
_Escala._aplicar = false;
double h[ordenNum + ordenDen + 1];
_respuestaImpulso(ordenNum + ordenDen + 1, h);
_Pade(h, num, ordenNum, den, ordenDen);
}
void FractionalOrderFilter::_Pade(double* h, double* num, int ordenNum, double* den, int ordenDen) {
double dH[ordenDen * ordenDen];
for (int i = 0; i < ordenDen; i++) {
for (int j = 0; j < ordenDen; j++) {
dH[ordenDen * i + j] = h[ordenNum + i - j];
}
}
for (int i = 0; i < ordenDen; i++) {
den[i] = -h[ordenNum + 1 + i];
}
_gaussJordan(dH,den,ordenDen);
std::copy(den, den + ordenDen, &den[1]);
den[0] = 1;
_numPade(h,num,ordenNum,den,ordenDen);
}
void FractionalOrderFilter::_numPade(double* h, double* num,double ordenNum,double* den,double ordenDen){
for (int n = 0; n <= ordenNum; n++) {
int limite = (n < ordenDen) ? n : ordenDen;
for (int k = 0; k <= limite; k++) {
num[n] += den[k] * h[n - k];
}
num[n] = num[n] * _Escala._valor;
}
}
void FractionalOrderFilter::_Shank(double* num, int ordenNum, double* den, int ordenDen, int cantidad) {
_ponerCeros(num,ordenNum);
_ponerCeros(den,ordenDen);
if (cantidad<=ordenNum+ordenDen){
return;
}else if(cantidad == ordenNum+ordenDen +1){
_Pade(num,ordenNum,den,ordenDen);
}else{
_Escala._aplicar = false;
double h[cantidad];
_respuestaImpulso(cantidad,h);
_Shank(h, num, ordenNum, den, ordenDen,cantidad);
_numPade(h, num, ordenNum, den, ordenDen);
}
}
void FractionalOrderFilter::_Shank(double* h, double* num, int ordenNum, double* den, int ordenDen, int cantidad) {
int sHTH = ordenDen * ordenDen;
double HTH[sHTH];
_ponerCeros(HTH,sHTH);
_MatrixParams params = {HTH, h, h, ordenDen, ordenNum, cantidad - 1, den, false};
_calcularMatriz(params);
_gaussJordan(HTH, den, ordenDen);
std::copy(den, den + ordenDen, &den[1]);
den[0] = 1.0;
}
void FractionalOrderFilter::_Prony(double* num, int ordenNum, double* den, int ordenDen, int cantidad) {
_ponerCeros(num,ordenNum);
_ponerCeros(den,ordenDen);
if (cantidad<=ordenNum+ordenDen){
return;
}else if(cantidad == ordenNum+ordenDen +1){
_Pade(num,ordenNum,den,ordenDen);
}else{
_Escala._aplicar = false;
double h[cantidad];
_respuestaImpulso(cantidad, h);
_Shank(h, num, ordenNum, den, ordenDen,cantidad);
_Prony(h, num, ordenNum, den, ordenDen,cantidad);
_numPade(h, num, ordenNum, den, -1);
}
}
void FractionalOrderFilter::_Prony(double* h, double* num, int ordenNum, double* den, int ordenDen, int cantidad) {
double g[cantidad];
g[0] = 1.0;
for (int n = 1; n < cantidad; n++) {
g[n] = 0.0;
for (int k = 0; k <= ((n < ordenDen) ? n : ordenDen); k++) {
g[n] -= den[k] * g[n - k];
}
}
int sGTG = (ordenNum + 1) * (ordenNum + 1);
double GTG[sGTG];
_ponerCeros(GTG,sGTG);
_MatrixParams params = {GTG, h, g, ordenNum + 1, 0, cantidad, num, true};
_calcularMatriz(params);
_gaussJordan(GTG, num, ordenNum + 1);
}
{0.0707106781186548,
0.0707106781186548,
0.0353553390593274,
0.0353553390593274,
0.0265165042944955,
0.0265165042944955,
0.0220970869120796,
0.0220970869120796,
0.0193349510480697,
0.0193349510480697,
0.0174014559432627,
0.0174014559432627,
0.0159513346146575,
0.0159513346146575,
0.0148119535707534,
0.0148119535707534,
0.0138862064725813,
0.0138862064725813,
0.0131147505574379,
0.0131147505574379,
0.012459013029566,
0.012459013029566,
0.0118926942554948,
0.0118926942554948,
0.0113971653281825,
0.0113971653281825,
0.0109588128155601,
0.0109588128155601,
0.0105674266435758,
0.0105674266435758,
0.01021517908879,
0.01021517908879,
0.00989595474226528,
0.00989595474226528,
0.00960489724984571,
0.00960489724984571,
0.00933809454846111,
0.00933809454846111,
0.00909235521823845,
0.00909235521823845,
0.00886504633778249,
0.00886504633778249,
0.00865397380593052,
0.00865397380593052,
0.00845729258306847,
0.00845729258306847,
0.00827343839648002,
0.00827343839648002,
0.00810107509655336,
0.00810107509655336,
0.00793905359462229,
0.00793905359462229,
0.0077863794870334,
0.0077863794870334,
0.00764218727431056,
0.00764218727431056,
0.00750571964441216,
0.00750571964441216,
0.00737631068502574,
0.00737631068502574,
0.00725337217360864,
0.00725337217360864,
0.00713638229984076,
0.00713638229984076,
0.00702487632640575,
0.00702487632640575,
0.00691843880630869,
0.00691843880630869,
0.0068166970591571,
0.00681669705915709,
0.00671931567259771,
0.00671931567259771,
0.00662599184381163,
0.00662599184381163,
0.00653645141348985,
0.00653645141348985,
0.00645044547383867,
0.00645044547383867,
0.0063677474549433,
0.0063677474549433,
0.00628815061175651,
0.00628815061175651,
0.0062114658481985,
0.0062114658481985,
0.00613751982619614,
0.00613751982619614,
0.00606615331658921,
0.00606615331658921,
0.00599721975617342,
0.00599721975617342,
0.00593058398110483,
0.00593058398110483,
0.00586612111174499,
0.00586612111174499,
0.00580371556800302,
0.00580371556800302,
0.00574326019750299,
0.00574326019750299,
0.0056846555016101,
0.0056846555016101}
{{0.0707106781186548,
0.0353553390600263,
-0.141421356236942,
-0.0618718433550409,
0.0928077650301804,
0.0331456303687676,
-0.0220970869118471,
-0.00552427172812658,
0.0013810679319814,
0.000138106793203062},{1,
-0.499999999990116,
-2.00000000000469,
0.874999999982521,
1.31250000000705,
-0.468749999990503,
-0.31250000000294,
0.0781249999983837,
0.0195312500002931,
-0.00195312499995776}}
{{0.0707106781186548,
0.00498890047350993,
-0.177026201633966,
-0.0104911888889076,
0.152112724513057,
0.00706521551752149,
-0.0506123254331939,
-0.00158673523758486,
0.00490139563555662,
6.03763413473617e-05},{1,
-0.929446292890327,
-2.07408225957115,
1.89043758989521,
1.3875254664944,
-1.22224353576288,
-0.318681155767505,
0.260419083129573,
0.0152884909573334,
-0.00918739826895705}}
{{0.0707106773187466,
0.00498607879616071,
-0.177024377758271,
-0.0104092897243263,
0.152098161778624,
0.00667325986292348,
-0.0505832096237445,
-0.000969870617454433,
0.00488490676159391,
-0.000247825369919612},{1,
-0.929446292890327,
-2.07408225957115,
1.89043758989521,
1.3875254664944,
-1.22224353576288,
-0.318681155767505,
0.260419083129573,
0.0152884909573334,
-0.00918739826895705}}
Compiling Program...
Standard Input:
Interactive Console
Text
Program is not being debugged. Click "Debug" button to start program in debug mode.
| # | Function | File:Line |
|---|
| Variable | Value |
|---|
| Register | Value |
|---|
| Expression | Value |
|---|
| # | Description |
|---|