question about dely lines

3 次查看（过去 30 天）
Muhsin Zerey
Muhsin Zerey 2024-12-1
评论： Muhsin Zerey 2025-1-14
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
  2 个评论
Paul
Paul 2024-12-25
编辑：Paul 2024-12-25
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
Muhsin Zerey 2025-1-13
Hi, yes they are all constants. I also already have the difference euqation correctly. Should be like this.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Walter Roberson
Walter Roberson 2024-12-1
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.
  8 个评论
显示 6更早的评论
Walter Roberson
Walter Roberson 2025-1-14
You need to initialize y(1) through y(dm0_prime)
Muhsin Zerey
Muhsin Zerey 2025-1-14
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)

请先登录，再进行评论。

Paul
Paul 2024-12-26
Ran in:
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>
  1 个评论
Muhsin Zerey
Muhsin Zerey 2025-1-13
Hi Paul, I reedited my question and hope you can understand my problem better.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Audio Processing Algorithm Design 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by