question about dely lines

EDITED:

Hey guys, I want to implement an allpass filter but i struggle with the difference equation and its implementation:

heres the structure

and here are the difference equations:

So finally I got the difference equation. I also tried to implemend it into my process function. (d(n) is a delay line in my code before i wanted to implemt the allpass, therefore I commented it out, but can be useful to compare). m(k) and m'(k) are both delays that are calculated. zeta is set to be one and is therefore not in the equation. The plugin sounds wrong and horrible if I try this way. Anyone got an Idea?

function out = process(plugin, in) 
            out = zeros(size(in));
            
            for i = 1:size(in,1)
                % Summieren der L/R - Kanle
                inL = in(i,1);
                inR = in(i,2);
                inSum = (inL + inR)/2;
                plugin.buffInput(plugin.pBuffInput + 1) = inSum;
               
               
                

                % loop over delay lines
                for n=1:plugin.N
                   % plugin.y_a = 0;
                    % d_n = gain * delayed v_n
                    for k=1:plugin.N
                        
                        
                        % if k == 2 && mod(plugin.pBuffDelayLines,2) == 0
                        %     plugin.gy(k) = 0;
                        % 
                        % end

                        plugin.Dg(k) = sqrt(1-plugin.g(k)^2);
                        
                        %plugin.d(k) = plugin.g(k)*plugin.buffDelayLines(k, mod(plugin.pBuffDelayLines + plugin.m(k), plugin.maxDelay +1) + 1); 
                        
                        % d(k) = (((sqrt(1-plugin.g(k)^2)^2)+ plugin.g(k)^2 + plugin.g(k)^2) * x1_m0p) + (plugin.g(k) * x1_m0) - (plugin.g(k) * y_m0p);
                        x1_m0p = plugin.buffDelayLines(k, mod(plugin.pBuffDelayLines + plugin.m(k)+plugin.m'(k)+1, plugin.maxDelay +1) + 1);
                        x1_m1p =plugin.buffDelayLines(k, mod(plugin.pBuffDelayLines+ plugin.m(k) +1, plugin.maxDelay +1) + 1);
                        plugin.d(k)= (plugin.Dg(k)^2+plugin.g(k)^2)*x1_m0p +  plugin.g(k)*x1_m1p- plugin.g(k)*plugin.y_a(k);
                        plugin.y_a(k) = plugin.d(k);
                    

                    end
                    %generate time variant matrix
                    %generateTIFDNmatrix(plugin,buffA);

                    % f_n = A(n,:) * d'
                    
                    plugin.f(n) = plugin.A(n,:) * plugin.d(:);
                    % v_n with pre delay
                    plugin.v(n) = plugin.b(n) * plugin.buffInput(mod(plugin.pBuffInput + plugin.preDelayS, (plugin.maxPreDelay * plugin.fs + 1)) + 1) ...
                        + plugin.f(n); %An pe delay noch arbeiten
                    
                   
                    

                    plugin.buffDelayLines(n, plugin.pBuffDelayLines + 1) = plugin.v(n);
                    % output lines
                   
                    plugin.s(n) = plugin.c(n)* plugin.d(n);
                    out(i,:) = out(i,:) + real(plugin.s(n));
                    
                    
                end

                % Assign to output
                out(i,1) = plugin.mix/100 * out(i,1) + (1.0 - plugin.mix/100) * in(i,1);
                out(i,2) = plugin.mix/100 * out(i,2) + (1.0 - plugin.mix/100) * in(i,2);

                calculatePointer(plugin);
            end 
end

[Paul]

Are g_0, D_g0, and zeta_0 all constants? Is delta^m0 an integer?

The input to the filter is x_1(n) and the output is y_1(n) ?

[Muhsin Zerey]

Hi, yes they are all constants. I also already have the difference euqation correctly. Should be like this.

Answer 1

[Walter Roberson]

you cannot implement those equations.

e(n) is defined in terms of d(n)

d(n) is defined in terms of e(n - something)

Substituting, e(n) is defined in terms of e(n - something)

This is infinite recursion, and so has no solution.

[Muhsin Zerey]

the structure above is essentially an allpass filter. Did I do the equations wrong? Is there another way you would do it maybe?

[Walter Roberson]

There are a few possibilities:

• the equations might be wrong
• you might be missing a termination condition for the recursion -- for example e(0) might be a specific value and potentially it might be provable that the recursive sequence always eventually leads to e(0)
• some or all of the given equations might be irrelevant. For example, the "real" d(n) might be a vector of values and the given equation form of d(n) might be irrelevant to the situation

[Muhsin Zerey]

Hi Walter, I edited this thread and would love to hear from you what you think.

[Walter Roberson]

@Muhsin Zerey

Unless is identical to zero, y[n] is defined using infinite recursion. For example if is 1, then y[0] is defined in terms of some stuff together with y[n-dm'] which would be y[0-1] which would be y[-1]. In turn y[-1] would be defined in terms of some stuff together with y[-1-1] which would be y[-2] . In turn y[-2] would be defined in terms of y[-3] and so on. The same kind of problem happens if dm' is negative, except then y[0] would be defined in terms of y[1] which would be defined in terms of y[2] and so on. You can only escape if dm' is 0, in which case y[n] would be defined in terms of some stuff together with g0*y[n] -- in which case you could isolate the y[n] term on the left as (1-g0)*y[n] = stuff leading to a direct definition of y[n] = stuff/(1-g0)

[Paul]

Typically one initializes the process with appropriate initial conditions and then iterates to update y[n].

[Muhsin Zerey]

I initiated y(n) = 0 from begin with. So the first time y(n) is 0.

[Walter Roberson]

You need to initialize y(1) through y(dm0_prime)

[Muhsin Zerey]

Everything is done with ring buffers. so the buffer delay line is a ring buffer and in the process function you save all the values of v(n) into the buffer delay line. So d(n) as well as y(n) is depented on v(n) and v(n) is pre delay + the outcome of the mixing matrix*d(n)

Answer 2

[Paul]

One can attack this symbolically if the parameters in the problems aren't known. If they are, one can proceed numerically using the Control System Toolbox. Example of the latter

Define the constants, assume a 4 sample delay

g_0    = 0.5;
D_g0   = sqrt(3)/2;
zeta_0 = 1;
delta_m0 = 4;

Define the lti objects for the three equations

sys1 = ss([g_0,D_g0*zeta_0],'Ts',-1,'InputDelay',[delta_m0,0],'InputName',{'x1','d'},'OutputName','y1');
sys2 = ss([zeta_0*D_g0,-g_0],'Ts',-1,'InputDelay',[delta_m0,0],'InputName',{'x1','d'},OutputName = 'e');
sys3 = ss(1,'Ts',-1,'InputDelay',delta_m0,'InputName','e','OutputName','d');

Connect all together

sys = connect(sys1,sys2,sys3,'x1',{'y1','e','d'});

With the selected constants, the system from x1 to y1 is allpass

opts = bodeoptions;
opts.MagUnits = 'abs';
bodeplot(sys(1,1),opts);

Plot the outputs with an input for x1

N = 50;
x1 = [ones(N/2,1);-ones(N/2,1)];
[z,k] = lsim(sys,x1);
y1 = z(:,1);e = z(:,2); d = z(:,3);
figure
hold on
stem(k,y1,'DisplayName','y1');
stem(k,e ,'DisplayName','e');
stem(k,d ,'DisplayName','d');
legend

Check that the outputs satisfy the original difference equations.

x1s = @(n) interp1(k,x1,n,'linear',0);
es  = @(n) interp1(k,e, n,'linear',0);
ds  = @(n) interp1(k,d, n,'linear',0);

[norm(y1 - ( g_0*x1s(k-delta_m0) + D_g0*zeta_0*ds(k) ));
norm(e - ( D_g0*zeta_0*x1s(k-delta_m0) - g_0*ds(k) ));
norm(d - es(k-delta_m0))]

ans = 3×1

1.0e-15 * 0.3168 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Muhsin Zerey]

Hi Paul, I reedited my question and hope you can understand my problem better.