To convert the sample rate filtering of dsp.FarrowRateConverter to a dsp.NewtonRateConverter by changing the Mathematical Implementation in the dsp.FarrowRateConverter ?

So my task is to override the sample rate filtering of system object dsp.FarrowRateConverter into dsp.NewtonRateConverter .

Theoretically Newton Structure which is the converted structure looks like this :

The mathematical representation of Farrow Rate structure :

Here it is the filter Transfer function is H

mu = [1 mu mu^2 mu^3]^T it is the delay vector at output

Z represents state vector [1 z^-1 z^-2 z^-3]^T

C is the matrix for Lagrange Polynomial

Basically through transformation matrices Tz,Td = Td1 cross Td2 the Filter transfer function changes into Newton rate Structure :

Td and Tz can be calculated as the following way :

so d^T is represented as the new delay vector for Newton Rate Converter

where

Aim is to make Farrow Rate Converter as Newton rate converter:

So I tried to override dsp.FarrowRateConversion I tried to do method overriding Here is my code:

classdef (StrictDefaults)NewtonRateConverter< dsp.FarrowRateConverter &...

dsp.internal.SupportScalarVarSize & ...

dsp.internal.MultirateEngine &...

dsp.internal.FilterAnalysisMultirate

%This new subclass inherits the properties of farrow rate converter

methods

function obj = NewtonRateConverter(varargin)

% Constructor for the NewtonRateConverter class

setProperties(obj, nargin, varargin{:}, ...

'InputSampleRate', ...

'OutputSampleRate', ...

'OutputRateTolerance', ...

'PolynomialOrder');

end

end

properties

TransformationDelayMatrix1

end

methods

function computeTransformationDelayMatrix1(obj)

% Access PolynomialOrder from superclass

M = obj.PolynomialOrder;

% Compute TransformationDelayMatrix1 using compute_Td1

obj.TransformationDelayMatrix1 = compute_Td1(M);

end

end

properties

TransformationDelayMatrix2

end

methods

function computeTransformationDelayMatrix2(obj)

% Access PolynomialOrder from superclass

M = obj.PolynomialOrder;

% Compute TransformationDelayMatrix2 using compute_Td2

obj.TransformationDelayMatrix2 = compute_Td2(M);

end

end

properties

TransformationStateDelayMatrix

end

methods

function computeTransformationStateDelayMatrix(obj)

% Access PolynomialOrder from superclass

M = obj.PolynomialOrder;

% Compute TransformationStateDelayMatrix using calculate_Tz

obj.TransformationStateDelayMatrix = calculate_Tz(M);

end

end

properties

CombinedTransformationDelayMatrix

end

methods

function computeCombinedTransformationDelayMatrix(obj)

% Ensure that both TransformationDelayMatrix1 and TransformationDelayMatrix2 are computed

if isempty(obj.TransformationDelayMatrix1)

obj.computeTransformationDelayMatrix1();

end

if isempty(obj.TransformationDelayMatrix2)

obj.computeTransformationDelayMatrix2();

end

% Compute CombinedTransformationDelayMatrix as the Kronecker product of Td1 and Td2

obj.CombinedTransformationDelayMatrix = kron(obj.TransformationDelayMatrix1, obj.TransformationDelayMatrix2);

end

end

properties

NewtonCoefficients

end

methods

function computeNewtonCoefficients(obj)

% Ensure all required matrices are computed

if isempty(obj.CombinedTransformationDelayMatrix)

obj.computeCombinedTransformationDelayMatrix();

end

if isempty(obj.TransformationStateDelayMatrix)

obj.computeTransformationStateDelayMatrix();

end

% Access PolynomialCoefficient (P) from superclass

P = obj.Coefficients;

% Compute (T_d^T)^-1 * P * T_z^-1

Td_transpose = obj.CombinedTransformationDelayMatrix';

Tz = obj.TransformationStateDelayMatrix;

% Use matrix division for better accuracy and performance

Td_inv_transpose_P = Td_transpose \ P; % (T_d^T)^-1 * P

obj.NewtonCoefficients = Td_inv_transpose_P / Tz; % Resulting matrix (T_d^T)^-1 * P * T_z^-1

end

% function computeNewtonCoefficients(obj)

% % Ensure all required matrices are computed

% if isempty(obj.CombinedTransformationDelayMatrix)

% obj.computeCombinedTransformationDelayMatrix();

% end

% if isempty(obj.TransformationStateDelayMatrix)

% obj.computeTransformationStateDelayMatrix();

% end

%

% % Access Polynomial Coefficient from superclass

% P = obj.Coefficients;

% Td_transpose = obj.CombinedTransformationDelayMatrix';

% if size(P, 1) ~= size(Td_transpose, 1)

% % If P is not in the right orientation, transpose it

% P = P';

% end

% % Compute the final transformation

% Td_inv_transpose_P = kron(pinv(Td_transpose ) , P);

% %final

% obj.NewtonCoefficients =Td_inv_transpose_P;

%

% end

end

properties

SplineCoefficients

end

methods

function computeSplineCoefficients(obj, u)

% Method to compute and store B-spline coefficients

% Inputs:

% s - Input sequence (1xN array)

% Outputs:

% obj.SplineCoefficients - Final B-spline coefficients (1xN array)

% Ensure input sequence length matches the polynomial order

M = obj.PolynomialOrder;

if length(u) ~= M

error('Input sequence length must match PolynomialOrder.');

end

% Calculate B-spline coefficients

obj.SplineCoefficients = calculateBsplineCoefficients(u);

end

end

properties

TransformedSplineCoefficients

end

methods

function computeTransformedSplineCoefficients(obj)

% Ensure necessary matrices are computed

if isempty(obj.TransformationStateDelayMatrix)

obj.computeTransformationStateDelayMatrix();

end

if isempty(obj.TransformationDelayMatrix1)

obj.computeTransformationDelayMatrix1();

end

% Access computed spline coefficients and transformation matrices

C_spline = obj.SplineCoefficients;

T_mu = obj.TransformationDelayMatrix1; % Assuming T_mu is the delay transformation matrix

T_z = obj.TransformationStateDelayMatrix;

% Apply the transformation: T_mu^(-T) * C_spline * T_z^(-1)

% Using matrix division instead of inv(A)*b for better performance and accuracy

T_mu_inv_transpose_C_spline = T_mu' \ C_spline; % Equivalent to inv(T_mu') * C_spline

obj.TransformedSplineCoefficients = T_mu_inv_transpose_C_spline / T_z; % Equivalent to T_mu^(-T) * C_spline * inv(T_z)

end

end

properties

ReconfigurableSplineCoefficients

end

methods

function computeReconfigurableSplineCoefficients(obj, cs)

% Method to compute ReconfigurableSplineCoefficients based on cs

% Inputs:

% cs - The weighting coefficient (scalar, in the range [0, 1])

if cs < 0 || cs > 1

error('The coefficient cs must be in the range [0, 1].');

end

% Calculate Q_Spl_cs as per the formula

obj.ReconfigurableSplineCoefficients = obj.Coefficients + cs * (obj.SplineCoefficients - obj.Coefficients);

end

end

properties

TransformedReconfigurableSplineCoefficients

end

methods

function computeTransformedReconfigurableSplineCoefficients(obj, cs)

% Method to compute TransformedReconfigurableSplineCoefficients

% by applying the transformation to ReconfigurableSplineCoefficients

% Ensure that ReconfigurableSplineCoefficients are computed

if isempty(obj.ReconfigurableSplineCoefficients)

obj.computeReconfigurableSplineCoefficients(cs);

end

% Ensure that Transformation matrices are computed

if isempty(obj.TransformationDelayMatrix1)

obj.computeTransformationDelayMatrix1();

end

if isempty(obj.TransformationStateDelayMatrix)

obj.computeTransformationStateDelayMatrix();

end

% Access the reconfigurable spline coefficients and transformation matrices

C_spline_reconfigurable = obj.ReconfigurableSplineCoefficients;

T_mu = obj.TransformationDelayMatrix1; % Assuming T_mu is the delay transformation matrix

T_z = obj.TransformationStateDelayMatrix;

% Apply the transformation: T_mu^(-T) * C_spline_reconfigurable * T_z^(-1)

% Using matrix division for efficiency and numerical stability

T_mu_inv_transpose = T_mu' \ eye(size(T_mu)); % Equivalent to inv(T_mu')

T_z_inv = T_z \ eye(size(T_z)); % Equivalent to inv(T_z)

obj.TransformedReconfigurableSplineCoefficients = T_mu_inv_transpose * C_spline_reconfigurable * T_z_inv;

end

end

% properties

% CoefficientType

% end

methods

function y = farrowloop(obj, u, coefficient)

% Newton-based rate conversion loop to replace the FarrowLoop

% Uses transformed delay, state, and coefficient matrices

if isempty(obj.NewtonCoefficients)

obj.computeNewtonCoefficients();

end

% Select coefficients based on CoefficientType

% switch obj.CoefficientType

switch coefficient

case "Newton"

if isempty(obj.NewtonCoefficients)

obj.computeNewtonCoefficients();

end

coefficients = obj.NewtonCoefficients;

case "Spline"

if isempty(obj.SplineCoefficients)

obj.computeSplineCoefficients(u);

end

coefficients = obj.SplineCoefficients;

case "TransformedSpline"

if isempty(obj.TransformedSplineCoefficients)

obj.computeTransformedSplineCoefficients();

end

coefficients = obj.TransformedSplineCoefficients;

case "ReconfigurableSpline"

if isempty(obj.ReconfigurableSplineCoefficients)

cs = 0.5; % Example value for cs, can be parameterized

obj.computeReconfigurableSplineCoefficients(cs);

end

coefficients = obj.ReconfigurableSplineCoefficients;

case "TransformedReconfigurableSplineCoefficients"

if isempty(obj.TransformedReconfigurableSplineCoefficients)

cs = 0.5; % Example value for cs, can be parameterized

obj.computeTransformedReconfigurableSplineCoefficients(cs);

end

coefficients = obj.TransformedReconfigurableSplineCoefficients;

otherwise

error('Invalid CoefficientType specified.');

end

% Transform delay vector (Tnext) and state vector (States)

Td = obj.CombinedTransformationDelayMatrix;

Tz = obj.TransformationStateDelayMatrix;

mu_transformed = (kron(Td , obj.Tnext))'; % (T_d * Tnext)^T

Z_transformed = kron(Tz , obj.States); % T_z * States

% Initialize output

y = zeros(size(u));

% Perform Newton-based interpolation using transformed coefficients

for k = 1:size(u, 1)

% Update transformed states with input and transformed coefficients

Z_transformed = Tz * [u(k, :); Z_transformed(1:end-1, :)];

obj.States = [u(k, :); obj.States(1:end-1, :)];

% Compute the polynomial interpolation using oefficients

intermediate_kron = kron(coefficients, Z_transformed);

y(k, :) = kron(mu_transformed', intermediate_kron);

%

% Update fractional delay and Tnext based on interpolation and decimation

obj.Tnext = obj.Tnext - obj.PrivM;

if obj.Tnext <= 0

obj.Tnext = obj.Tnext + obj.PrivL;

end

end

end

end

% methods(Access=protected)

% function [y, lTnext, Z] = farrowLoop(obj,fd,MultiplicandPrototype, u, y)

% nu = size(u,1); % Input frame size

% nc = size(u,2); % Number of channels

%

% L = coder.const(obj.PrivL);

% M = coder.const(obj.PrivM);

% lTnext = obj.Tnext; % Phase

% C_dbl = coder.const(obj.PolynomialCoefficients);

% Np = coder.const(size(C_dbl,1) - 1);

% C = coder.const(obj.PrivCoeff);

% Z = obj.States;

% Linv = coder.const(obj.PrivLinv);

% yIdx = int32(1);

%

% % Set up multiplicand

% pMpc = zeros([1 nc], 'like', MultiplicandPrototype);

%

% for k = 1:nu

% % Update internal states and apply filter coefficients

% Z(:) = Tz *[u(k,:); Z(1:end-1,:)];

%

% % Obtain polynomial coefficients for the current interval

% Cz = C*Z;

%

% % Loop fractional delay over current interval

% while (lTnext >0)

% % Compute fractional delay

% fd(:) = lTnext .* Linv;

%

% % Apply polynomial using Horner's rule

% pMpc(:) = Cz(1,:);

% yAcc = fd * pMpc;

%

% for pIdx = 2:Np

% pMpc(:) = Cz(pIdx,:) + yAcc;

% yAcc = fd * pMpc;

% end

% y(yIdx, :) = yAcc + Cz(Np+1,:);

%

% % Update local loop states.

% lTnext(:) = lTnext - M;

% yIdx(:) = yIdx + 1;

% end

%

% lTnext(:) = lTnext + L;

% end

% end

% function [y, lTnext, Z] = farrowLoop(obj, fd, ~, u, y, coefficientType)

% % Generalized loop for Newton and Spline-based interpolation

% nu = size(u, 1); % Input frame size

% nc = size(u, 2); % Number of channels

%

% % Initialize constants and variables

% L = coder.const(obj.PrivL);

% M = coder.const(obj.PrivM);

% lTnext = obj.Tnext; % Phase

% Z = obj.States;

% Linv = coder.const(obj.PrivLinv);

% yIdx = int32(1);

%

% % Select the coefficient set based on input type

% switch coefficientType

% case "Newton"

% coefficients = obj.NewtonCoefficients;

% case "Spline"

% coefficients = obj.SplineCoefficients;

% case "TransformedSpline"

% coefficients = obj.TransformedSplineCoefficients;

% case "ReconfigurableSpline"

% coefficients = obj.ReconfigurableSplineCoefficients;

% otherwise

% error('Invalid coefficient type specified.');

% end

%

% % Precalculate transformations

% Td_transpose = obj.CombinedTransformationDelayMatrix';

% Tz = obj.TransformationStateDelayMatrix;

%

% % Apply the transformations to the coefficients

% transformedCoefficients = Td_transpose \ coefficients / Tz;

%

% % Loop for interpolation

% for k = 1:nu

% % Update states

% Z = [u(k, :); Z(1:end-1, :)];

%

% % Get the polynomial coefficients for this interval

% Cz = transformedCoefficients * Z;

%

% % Perform the filtering step using the transformed coefficients

% while lTnext > 0

% % Compute fractional delay

% fd(:) = lTnext .* Linv;

%

% % Perform interpolation using Horner's rule on transformed coefficients

% yAcc = Cz(1, :);

% for pIdx = 2:size(Cz, 1)

% yAcc = fd .* yAcc + Cz(pIdx, :);

% end

%

% % Store result

% y(yIdx, :) = yAcc;

% yIdx = yIdx + 1;

%

% % Update phase

% lTnext = lTnext - M;

% end

%

% % Increment for next phase

% lTnext = lTnext + L;

% end

% end

% end

methods(Static)

function c = calculateBsplineCoefficients(s)

% Function to calculate B-spline coefficients for a given sequence s

% Inputs:

% u - Input sequence to be interpolated (1xN array)

% Outputs:

% c - Final B-spline coefficients (1xN array)

% Parameters

N = length(u); % Length of input sequence

iterations = 10; % Number of iterations for alpha convergence

alpha_values = zeros(1, iterations); % Initialize array to store alpha values

% Initial condition for alpha

alpha_values(1) = -1/4;

% Iteratively calculate alpha values

for i = 2:iterations

alpha_values(i) = -1 / (4 + alpha_values(i - 1));

end

% Use the last value as the converged alpha

alpha = alpha_values(end);

% Value for cubic B-splines

J = min(10, N);

b_i = -alpha / (1 - alpha^2);

% Initialize c_plus array

c_plus = zeros(1, N);

% Step 1: Calculate c^+(1) using Equation (2)

c_plus(1) = sum(alpha.^(0:J-1) .* s(1:J));

% Step 2: Calculate c^+(k) for k = 2, ..., N using Equation (3)

for k = 2:N

c_plus(k) = s(k) + alpha * c_plus(k - 1);

end

% Calculate c^-(N) using Equation (4)

c_minus = zeros(1, N);

c_minus(N) = b_i * (2 * c_plus(N) - s(N));

% Calculate c^-(k) for k = N-1, ..., 1 using Equation (5)

for k = N-1:-1:1

c_minus(k) = alpha * (c_minus(k + 1) - c_plus(k));

end

% Calculate final B-spline coefficients c(k) using Equation (6)

c = 6 * c_minus;

end

function Td1 = compute_Td1(M)

% Create a matrix of binomial coefficients

[I, J] = ndgrid(1:M, 1:M);

binomials = zeros(M, M);

for i = 1:M

binomials(i, 1:i) = arrayfun(@(ii, jj) nchoosek(ii-1, jj-1), I(i, 1:i), J(i, 1:i));

end

% Create a matrix for powers

powers = ((- (M - 1) / 2) .^ (I - J)) .* (J <= I);

% Element-wise multiplication to get Td1

Td1 = binomials .* powers;

end

function Td2 = compute_Td2(M)

% Initialize Td2 with identity

Td2 = eye(M + 1);

for n = 2:M+1

for k = 2:n-1

Td2(n, k) = Td2(n-1, k-1) - (n-2) * Td2(n-1, k);

end

end

Td2 = Td2(2:M+1, 2:M+1);

end

function T_z = calculate_Tz(M)

% Initialize the matrix T_z with zeros

T_z = zeros(M, M);

% Fill the matrix T_z using the given formula

for i = 1:M

for j = 1:i % j goes from 1 to i

T_z(i, j) = nchoosek(i-1, j-1) * (-1)^(j+1);

end

end

end

end

end

This is not getting overridden please could you give suggestions or debug the ccode ??? I tried to ovveride the farrow structure as a Newton structure and then I tried to use Spline Coefficients,Transformed Spline Coefficients ,Reconfigurable spline Coefficients ,Newton Coefficients

Then calculate execute the frequency response using fvtool see the cost and all parameterss.