# include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
/* Implements an IIR-Filter with rational transfer function,
whose nominator coeff. in A_in, denominator coeff. in B_in,
and input samples in x_in. This can be samples of a signal with M channels.
´ y_out is the filter output. It is assumed that A_in and B_in have the same
number of coefficients and that B_in[0]=1.
Usage: y=myfilter(A,B,x);
Remarks: 1) Use mex with the C++ compiler to compile the routine
2) This filter function is a so-called "Direct Form II Transposed"
implementation
3) Nothing is done to correct the phase distortions that come
along with this IIR-filter. This is left for the user's
own work.
R. Brigola, TH Nürnberg, December 2014
*/
#define A_in prhs[0]
#define B_in prhs[1]
#define x_in prhs[2]
#define y_out plhs[0]
double *A, *B, *x, *y;
int n, m, N, M;
double s;
int c , p, k;
/* The inputs must be noncomplex scalar double */
n = mxGetNumberOfElements(A_in); /* number of nominator coeff. of the filter */
m = mxGetNumberOfElements(B_in); /* number of denominator coeff. of the filter */
if(n!=m) /* Check the number of filter coefficients */
mexErrMsgTxt("ERROR: Unequal number of filter coefficients in nominator and denominator.");
else
N = mxGetM(prhs[ 2 ]); /* number of samples to be filtered */
M = mxGetN(prhs[ 2 ]); /* number of input channels */
/* Create matrix for the return argument */
y_out = mxCreateDoubleMatrix (N , M , mxREAL );
/* Assign pointers to input and output */
A = mxGetPr (A_in);
B = mxGetPr (B_in);
x = mxGetPr (x_in);
y = mxGetPr (y_out);
if(B[0]!=1) /* Check the normalization of B */
mexErrMsgTxt("ERROR: Filter coefficient B[0] is not 1.");
else
/* Filter algorithm, i.e. implement the corresponding difference equation */
/* Observe, that a (MxN)-matrix element M(m+1,n+1) here is addressed as M[n+N*m], */
/* when 0<=m<=M-1, 0<=n<=N-1 */
for (c= 0; c