- You will see updates in your activity feed.
- You may receive emails, depending on your notification preferences.

125 views (last 30 days)

Show older comments

Hello everybody!

How to find the coefficients of this difference equation in matlab in order to filter a music signal?

Thanks in advance!

Mathieu NOE
on 23 Feb 2021

helo

this looks like an echo flter ... but the answer is in the question , you have written all coefficients of the equation - !

pretty much what is done here :

infile='DirectGuitar.wav';

outfile='out_echo.wav';

% read the sample waveform

[x,Fs] = audioread(infile);

% normalize x to +/- 1 amplitude

x = x ./ (max(abs(x)));

% parameters

N_delay=20; % delay in samples

C = 0.7; % amplitude of direct sound

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y = zeros(length(x),1); % create empty out vector

y(1:N_delay)=x(1:N_delay); % to avoid referencing of negative samples

% for each sample > N_delay

for i = (N_delay+1):length(x)

y(i) = C*x(i) + (1-C)*(x(i-N_delay)); % add delayed sample

end

% write output

% normalize y to +/- 1 amplitude

y = y ./ (max(abs(y)));

audiowrite(outfile, y, Fs);

figure(1)

hold on

plot(x,'r');

plot(y,'b');

title('Echoed and original Signal');

sound(y,Fs);

Mathieu NOE
on 23 Feb 2021

the delay in samples = round( Fs x 0.15 sec)

do you know Fs ?

Mathieu NOE
on 24 Feb 2021

here is it : 3 filters in series :

infile='DirectGuitar.wav';

outfile='out_reverb.wav';

% read the sample waveform

[x,Fs] = audioread(infile);

% normalize x to +/- 1 amplitude

x = x ./ (max(abs(x)));

% parameters (3 delay filters in series)

N1_delay=15; % delay in samples

C1 = 0.7; % amplitude of direct sound

N2_delay=20; % delay in samples

C2 = 0.6; % amplitude of direct sound

N3_delay=50; % delay in samples

C3 = 0.5; % amplitude of direct sound

N_delay = max([N1_delay N2_delay N3_delay]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initialization %

y1 = zeros(length(x),1); % create empty out vector

y2 = y1;

y3 = y1;

y1(1:N_delay)=x(1:N_delay); % to avoid referencing of negative samples

y2(1:N_delay)=x(1:N_delay); % to avoid referencing of negative samples

y3(1:N_delay)=x(1:N_delay); % to avoid referencing of negative samples

% for each sample > N_delay

for i = (N_delay+1):length(x)

y1(i) = C1*x(i) + (1-C1)*(x(i-N1_delay)); % add delayed sample

y2(i) = C2*y1(i) + (1-C2)*(y1(i-N2_delay)); % add delayed sample

y3(i) = C3*y2(i) + (1-C3)*(y2(i-N3_delay)); % add delayed sample

end

% write output

% normalize y to +/- 1 amplitude

y3 = y3 ./ (max(abs(y3)));

audiowrite(outfile, y3, Fs);

figure(1)

hold on

plot(x,'r');

plot(y3,'b');

title('Echoed and original Signal');

sound(y3,Fs);

Mathieu NOE
on 24 Feb 2021

have you tried to do a fft of it ?

example below :

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load signal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% data

[data,Fs] = audioread('test_voice.wav');

channel = 1;

signal = data(:,channel);

samples = length(signal);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% FFT parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NFFT = 4096; %

OVERLAP = 0.75;

% spectrogram dB scale

spectrogram_dB_scale = 80; % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% options

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if you are dealing with acoustics, you may wish to have A weighted

% spectrums

% option_w = 0 : linear spectrum (no weighting dB (L) )

% option_w = 1 : A weighted spectrum (dB (A) )

option_w = 0;

%% decimate (if needed)

% NB : decim = 1 will do nothing (output = input)

decim = 4;

if decim>1

signal = decimate(signal,decim);

Fs = Fs/decim;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 1 : averaged FFT spectrum

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

% convert to dB scale (ref = 1)

sensor_spectrum_dB = 20*log10(sensor_spectrum);

% apply A weigthing if needed

if option_w == 1

pondA_dB = pondA_function(freq);

sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;

my_ylabel = ('Amplitude (dB (A))');

else

my_ylabel = ('Amplitude (dB (L))');

end

figure(1),plot(freq,sensor_spectrum_dB,'b');grid

title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);

xlabel('Frequency (Hz)');ylabel(my_ylabel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 2 : time / frequency analysis : spectrogram demo

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));

% FFT normalisation and conversion amplitude from linear to dB (peak)

sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg)); % NB : X=fft(x.*hanning(N))*4/N; % hanning only

% apply A weigthing if needed

if option_w == 1

pondA_dB = pondA_function(fsg);

sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));

my_title = ('Spectrogram (dB (A))');

else

my_title = ('Spectrogram (dB (L))');

end

% saturation of the dB range :

% saturation_dB = 60; % dB range scale (means , the lowest displayed level is XX dB below the max level)

min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;

sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

% plots spectrogram

figure(2);

imagesc(tsg,fsg,sg_dBpeak);colormap('jet');

axis('xy');colorbar('vert');grid

title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);

xlabel('Time (s)');ylabel('Frequency (Hz)');

function pondA_dB = pondA_function(f)

% dB (A) weighting curve

n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));

r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));

pondA = n./r;

pondA_dB = 20*log10(pondA(:));

end

function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)

% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).

% Linear averaging

% signal - input signal,

% Fs - Sampling frequency (Hz).

% nfft - FFT window size

% Overlap - buffer overlap % (between 0 and 0.95)

signal = signal(:);

samples = length(signal);

% fill signal with zeros if its length is lower than nfft

if samples<nfft

s_tmp = zeros(nfft,1);

s_tmp((1:samples)) = signal;

signal = s_tmp;

samples = nfft;

end

% window : hanning

window = hanning(nfft);

window = window(:);

% compute fft with overlap

offset = fix((1-Overlap)*nfft);

spectnum = 1+ fix((samples-nfft)/offset); % Number of windows

% % for info is equivalent to :

% noverlap = Overlap*nfft;

% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows

% main loop

fft_spectrum = 0;

for i=1:spectnum

start = (i-1)*offset;

sw = signal((1+start):(start+nfft)).*window;

fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only

end

fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum % Select first half

if rem(nfft,2) % nfft odd

select = (1:(nfft+1)/2)';

else

select = (1:nfft/2+1)';

end

fft_spectrum = fft_spectrum(select);

freq_vector = (select - 1)*Fs/nfft;

end

Mathieu NOE
on 25 Feb 2021

this will do averaged fft and spectrogram analysis of any signal , and in case it has some periodic content, it will show up as a peak in the spectrum (at the dominant frequency)

otherwise , if you are sure that your signal is basically a decaying (damped) sinus wave, you can also simply compute the time intervals (= period of oscillation) between consecutives crossing points (positive or negative slope of signal) , given a certain threshold , see example on steady sinus wave attached

RoBoTBoY
on 25 Feb 2021

Hello again!

I have these coefficients from a filter for dereverberation:

b = [6 0 0 0 0 0 0];

a = [1 0 2.45 0 2 0 0.55];

and this is the difference equation:

dereverb(n) = 6*reverb(n)-2.45*dereverb(n-2)-2*dereverb(n-4)-0.55*dereverb(n-6);

also I have these coefficients from a filter for reverb:

b = [0.1664 0 0.4084 0 0.3341 0 0.0911];

a = [1 0 0 0 0 0 0];

and this is the difference equation:

y(n) = 0.1664*x(n)+0.4084*x(n-2)+0.3341*x(n-4)+0.0911*x(n-6);

Have I calculated the dereverberation coefficients correctly? If yes, how to filter a reverb signal through dereverberation filter?

Thanks

Mathieu NOE
on 26 Feb 2021

hello

IMO, your dereverb filter is incorrect. I understand that you take the reverb filter (a FIR filter) and inverse numerator / denominator to create the dereverb filter. But that does not generate a causal , stable and realizable dereverb filter.

It is unstable as you can see by generating the impulse response :

b = [6 0 0 0 0 0 0];

a = [1 0 2.45 0 2 0 0.55];

dimpulse(b,a)

Mathieu NOE
on 26 Feb 2021

the usual technique used in acoustic room compensations , is to make a time delayed version of the flipped FIR filter

you can do that by creating the flipped FIR :

b_flipped = flipud(b')'

then you will see that putting the two filters in series lead to a perfectly flat tranfer function, that is , the impulse response of both filters in series is a sharp pulse (Dirac)

plot(conv(flipud(b')',b))

for your dereverb simulation, you only need to rewritte the difference equation based on b_flipped coefficients (as for the reverb case)

RoBoTBoY
on 28 Feb 2021

I wrote this:

b_dereverb = [4.6297 0 0 0 0 0 0];

a_dereverb = [1 0 2 0 1.3333 0 4];

flipped_b_dereverb = flipud(b_dereverb')';

flipped_a_dereverb = flipud(a_dereverb')';

y_derevereration = filter(flipped_b_dereverb,flipped_a_dereverb,reverb_piano);

y_derevereration = y_derevereration./(max(abs(y_derevereration))); % normalization of dereverb piano

figure(41);

plot(y_derevereration);

title('Dereverb Signal');

xlabel('Time');

pause(3);

sound(y_derevereration,Fs_piano);

But I don't have the expectantly results. I got the same reverb signal.

Why that?

Mathieu NOE
on 1 Mar 2021

hello

the "flip" technique applies only on FIR filter , not IIR; you can see that in many applications dealing with room acoustcs compesntion (electronic equalization of loudspeakers)

infile='DirectGuitar.wav';

outfile1='out_reverb2.wav';

outfile2='out_reverb_derev.wav';

% read the sample waveform

[x,Fs] = audioread(infile);

% normalize x to +/- 1 amplitude

x = x ./ (max(abs(x)));

%%%%%% reverb using FIR filter %%%%%%%%%%%%%%

% b = [0.1664 0 0.4084 0 0.3341 0 0.0911];

% a = [1 0 0 0 0 0 0];

% and this is the difference equation:

y = x;

for n = 7:length(x)

y(n) = 0.1664*x(n)+0.4084*x(n-2)+0.3341*x(n-4)+0.0911*x(n-6);

end

% write output

% normalize y to +/- 1 amplitude

y = y ./ (max(abs(y)));

audiowrite(outfile1, y, Fs);

%%%%%% dereverb using flipped FIR filter %%%%%%%%%%%%%%

% b_flipped = flipud(b')'

% b = [0.0911 0 0.3341 0 0.4084 0 0.1664];

% a = [1 0 0 0 0 0 0];

% and this is the difference equation:

y2 = y;

for n = 7:length(y)

y2(n) = 0.0911*y(n)+0.3341*y(n-2)+0.4084*y(n-4)+0.1664*y(n-6);

end

% write output

% normalize y to +/- 1 amplitude

y2 = y2 ./ (max(abs(y2)));

audiowrite(outfile2, y2, Fs);

figure(1)

hold on

plot(x,'r');

plot(y2,'b');

title('Echoed and original Signal');

sound(y2,Fs);

Mathieu NOE
on 2 Mar 2021

should be working, I tested it on another wav file and there was some improvements, even though the reverb + dereverb process was creating some masking effects - the final wav file sound is not as clean as the original one.

maybe there are more advanced methods...as attached papers

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

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

Start Hunting!Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)