Trying a wiener noise cancellation code but not able to filter out the noise

Question

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)

Answer 1

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