1
$\begingroup$

I am performing noise cancellation using Wiener filter for a project. The project statement is "to record our voice and take it as desired signal, add a noise, use wiener filter to estimate the previously added noise and finally plot the mean square error and estimated signal". So far, the algorithm looks correct but i am still not able to extract the filtered desired signal properly. The final output is the same as input(desired signal+noise).

clc;  
close all;
%taking the recorded signal as d
[d, fs] = audioread('ADSP recording .m4a'); 
%taking noisy signal 2 (v2 in Monson hayes diagram) as g 
n=length(d);
v2=0.5*randn(n,2);
p=20;
v2filt=filter(1, [1 -0.5] ,v2);
%%randn(1,n) is noisy signal 1 (v1) as per monson hayes diagram. we want to
% to approximate v1 using v2
v1=0.1*randn(n,2);
x=d+v1;
Rv1=covar(v2filt,p) ;
rxv=convm(x,p)'*convm(v2filt,p)/(n-1);
w=rxv(1,:)/Rv1;
v1hat=filter(w,1,v2filt);
appx=x-v1hat;
error=mse(d-appx);
%frequency domain mse
P = periodogram(d-appx,[],[],fs);
Hmss = dspdata.msspectrum(P,'Fs',fs,'spectrumtype','onesided'); 
%display
subplot(5,1,1),plot(d);
title('desired signal:');
subplot(5,1,2),plot(x);
title('desired + noise signal:');
subplot(5,1,3),plot(v1hat);
title('estimated signal:');
subplot(5,1,4),plot(appx);
title('Noise-removed signal');
subplot(5,1,5),plot(Hmss);
%title('frequency domain mean square error signal');

d = desired signal, v1 = added noise, v2 = secondary noise which is to be approximated to v1 using wiener filter, appx = the final cleaned signal (should contain only my desired signal)

$\endgroup$
1
  • $\begingroup$ command filter not used correctly: v2filt=filter(..) expects 1st input field to be a system, not a signal. The transfer function of the filter should be the system to input in 1st field of command filter. Obtain the transfer function with tf(1,[1 -.5]) $\endgroup$ Dec 15, 2022 at 20:42

1 Answer 1

1
$\begingroup$

I got a satisfactory output on the code by rectifying these errors:

  1. If you are recording a speech signal with your phone, make sure you know how many channels the recorded signal has. Mine was a dual channel recorder (Most of us dont know our phone recorder specs, but for stereo effect, dual channel recording is done these days by most mobiles). if your signal x(n) is a dual channel, that means its matrix representation has 2 columns. If you want a single channel, average both of them and store them in another array. Since both signals are your speech, it shouldn't be a problem.

  2. Nextly, in the above code, i have taken two different noise sources. When we take it that way, we are indirectly making them uncorrelated and we would never get a non-trivial solution to Wiener Hopf equation. I made a big mistake by taking those two uncorrelated noise sources. Instead, in my corrected version, I have taken a single randn() function as noise. For controllability of filter error, i have modelled two filter() functions v1 and v2 as AR(p) processes. That solved my problem instantly.

Here is the corrected Code ("ADSP recording.wav is my speech signal") . Also, Covar.m and convm.m are two matlab functions i created to create a covariance matrix and a convolutional matrix (These are also available in "Statistical Signal Processing" by Monson Hayes :

Main Code:

    clc;
clear all;
close all;
%taking the recorded signal as d
[d, fs] = audioread('ADSP recording.wav');
d=(d(:,1)+d(:,2))/2;
%taking noisy signal 2 (v2 in Monson hayes diagram) as g 
n=length(d);
g=randn(n,1);
v2=filter(1, [1 -0.3 0.5] ,g);
v1=filter(1, [1 -0.5 0.1] ,g);
%%randn(1,n) is noisy signal 1 (v1) as per monson hayes diagram. we want to
% to approximate v1 using v2
x=d+v1;
%sound(x,fs)
Rv2=covar(v2,5) ;
rxv=convm(x,5)'*convm(v2,5)/(n-1);
ww=rxv(1,:)/Rv2;
v1hat=filter(ww,1,v2);
appx=x-v1hat;
error=mse(d-appx)
sound(appx,fs)

As for the two defined functions conv.m and covar.m, here are their codes:

Conv.m

function X = convm(x,p)

N=length(x)+2*p-2;
x=x(:);
xpad = [zeros(p-1,1);x;zeros(p-1,1)];
for i=1:p
    X(:,i)=xpad(p-i+1:N-i+1);
end
end

covar.m

function R = covar(x,p)

x=x(:);
m=length(x);
x=x-ones(m,1)*(sum(x)/m);
R=convm(x,p).'*(convm(x,p)/(m-1));

end
$\endgroup$
1
  • $\begingroup$ Thanks for revisiting and posting your answer. $\endgroup$
    – Peter K.
    Dec 16, 2022 at 12:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.