How do I implement a low pass filter difference equation to filter sample by sample, samplewise
显示 更早的评论
Hi Guys, Im currently working on a reverberation alogrithm.
Heres the structure:
My supervisor told me that I have to implement a difference equation to filter the samples samplewise by just putting the samples in the equation and feed the result back to the system etc...
Now my question is how and where in the loop I can do that? The Loop Filters in the picture should be all second order low pass filters. How do I get the difference equation for such a filter?
Heres my code:
in = [ 1; 0 ]; % Dirac Impulse
Fs = 44100;
in = [in; zeros(3*Fs,1)]; % Space for reverb
% Define delay parameters
%Hadamard Matrix as feedback matrix
H = 0.75*hadamard(16);
delayTime = [0.02, 0.03, 0.05, 0.07];% in seconds 0.01 = 10ms 0.02 = 20ms 0.03 = ....
%feedbackGain = [0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65]; % adjust to control feedback level 0.5,0.53,0.55,0.58,0.6,0.63,0.65,0.68,0.7,0.72,0.75,0.77,0.8,0.81,0.85,0.88
feedbackGain =[0.1,-0.15,-0.2,0.25,-0.3,0.35,0.4,-0.45,-0.5,0.55,0.6,-0.65,-0.7,-0.75,-0.8,0.85]; %0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85
wetGain = [0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7]; % adjust to control wet signal level
% Convert delay time to samples
delaySamples = [101,233, 439,613, 853,1163, 1453,1667 , 1801,2053, 2269,2647, 3001,3607, 4153,4591]; %0.02s, 0.03, 0.04, 0.05 881, 1009, 1321,1511, 1777, 1999, 2203, 2351
%101, 439, 853, 1453 , 1801, 2269, 3001, 4153
%Summe delays = 30644
% delaySamples = [1200, 1800, 2200,2500];
% Create an empty array to store delayed audio
delayedAudio = zeros(size(in));
delayedAudio2 = zeros(size(in));
delayedAudio3 = zeros(size(in));
delayedAudio4 = zeros(size(in));
delayedAudio5 = zeros(size(in));
delayedAudio6 = zeros(size(in));
delayedAudio7 = zeros(size(in));
delayedAudio8 = zeros(size(in));
delayedAudio9 = zeros(size(in));
delayedAudio10 = zeros(size(in));
delayedAudio11 = zeros(size(in));
delayedAudio12 = zeros(size(in));
delayedAudio13 = zeros(size(in));
delayedAudio14 = zeros(size(in));
delayedAudio15 = zeros(size(in));
delayedAudio16 = zeros(size(in));
delayedAudios =[delayedAudio,delayedAudio2,delayedAudio3,delayedAudio4,delayedAudio5,delayedAudio6,delayedAudio7,delayedAudio8,delayedAudio9,delayedAudio10,delayedAudio11,delayedAudio12,delayedAudio13,delayedAudio14,delayedAudio15,delayedAudio16];
delayedAudios(1,:) = 1;
d = designfilt('lowpassfir','FilterOrder',2 ,'CutoffFrequency',2000,'SampleRate', Fs);
% Apply delay effect
for i = 1:length(in)-2*max(delaySamples)
for k = 1:16
if delayedAudios(i+delaySamples(k),k) == 0
delayedAudios(i+delaySamples(k),k) = filterfeedbackGain(k)*delayedAudios(i,k);
end
end
for k = 1:16
for j = 1:16
if delayedAudios(i+delaySamples(j)+delaySamples(k),k) == 0 && k~= j
delayedAudios(i+delaySamples(j)+delaySamples(k),k) = H(j,k)*delayedAudios(i+delaySamples(j),j);
end
end
end
end
%Original and Mix together
for k = 1:16
delayedAudios(:,k) = wetGain(k)*delayedAudios(:,k);
end
%outputAudioM = delayedAudios(:,16) + delayedAudios(:,15) + delayedAudios(:,14) + delayedAudios(:,13) + delayedAudios(:,12) +delayedAudios(:,11) +delayedAudios(:,10) +delayedAudios(:,9)+delayedAudios(:,8) +delayedAudios(:,7) +delayedAudios(:,6) +delayedAudios(:,5) +delayedAudios(:,4) +delayedAudios(:,3) +delayedAudios(:,2 )+delayedAudios(:,1);
outputAudioM = sum(delayedAudios,2);
outputAudioM(1) = outputAudioM(1)*0.02;
% %Write the output audio to a new file and plot
audiowrite('output_audio_with_delay.wav', outputAudioM, Fs);
plot(outputAudioM);
采纳的回答
hello
see this second order recursion - here used with a step input (works for any input signal)
I used the b,a coefficients of a second order low pass Butterworth filter (up to you to change filter type and cut off frequency)
all this can by simply implemented using gain and delay blocks in simulink
[b,a] = butter(2,0.1);
b0 = b(1);
b1 = b(2);
b2 = b(3);
a1 = a(2);
a2 = a(3);
% step response (step input at sample = 11)
% manual for loop coding IIR filter
x = [zeros(10,1); ones(60,1)];
samples = length(x);
y(1) = b0*x(1) + 0 + 0 + 0 + 0;
y(2) = b0*x(2) + b1*x(1) + 0 - a1*y(1) - 0;
for k = 3:samples
y(k) = b0*x(k) + b1*x(k-1) + b2*x(k-2) - a1*y(k-1) - a2*y(k-2);
end
figure(1)
plot(x,'k')
hold on
plot(y,'r')
6 个评论
Muhsin Zerey
2024-3-1
编辑：Muhsin Zerey
2024-3-1
Mathieu NOE
2024-3-1
编辑：Mathieu NOE
2024-3-1
hello
simply because matlab cannot work with zero or negative index , so if you have a recursion like above
y(k) = b0*x(k) + b1*x(k-1) + b2*x(k-2) - a1*y(k-1) - a2*y(k-2);
you must take care of how to generate the first 2 outputs (as for k = 1 and 2 , k - 2 is <= 0 )
Muhsin Zerey
2024-3-1
Mathieu NOE
2024-3-1
the recursion is valid for any time delay , but if you keep the same a and b coefficients , this would correspond to a different cut off frequency of your digital filter
if you need to have the same exact filtering for different sampling rates , then you have to compute a and b coefficients for each sampling rate and use them accordingly on the corresponding data
Muhsin Zerey
2024-3-5
the same logic, with minor modifications, can be used for a multichannel x input
here I generated 16 steps , with some time delta between the channels
clc
clearvars
[b,a] = butter(2,0.1);
b0 = b(1);
b1 = b(2);
b2 = b(3);
a1 = a(2);
a2 = a(3);
% step response (step input at sample = 11)
% manual for loop coding IIR filter
% multi channel input (70 samples / 16 channels
channels = 16;
for ci = 1:channels
x(:,ci) = [zeros(10+ci,1); ones(60-ci,1)];
end
% filtering the multi channel input
samples = size(x,1);
y(1,:) = b0*x(1,:) + 0 + 0 + 0 + 0;
y(2,:) = b0*x(2,:) + b1*x(1,:) + 0 - a1*y(1,:) - 0;
for k = 3:samples
y(k,:) = b0*x(k,:) + b1*x(k-1,:) + b2*x(k-2,:) - a1*y(k-1,:) - a2*y(k-2,:);
end
figure(1)
plot(x,'k')
hold on
plot(y,'r')
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Get Started with DSP System Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
