信号处理（四）Interference cancellation 消除干扰

You should download the following signals:
x1.wav and x2.wav

We consider two signals 𝑥1(𝑛) and 𝑥2(𝑛). We assume that they were constructed as follows:

x(n) is the useful signal that you should recover, w(n) is an additive interfering signal uncorrelated with x(n). The coefficients, a1, a2,b1 and b2 are real positive numbers aiming to create, on each signal, an unfavorable signal to interference ratio, that prevents you from listening to the original signal.

In order to recover that original signal, you should implement a dedicated power inversion technique. For this, you will linearly combine the two signals using a coefficient α through a linear equation as following:

You should determine the value of α that minimizes the power of the signal ε(n).

You will try to demonstrate that the signal to interference ratio, is reversed on ε(n), when compared to that of the signal x2(n). For this you must demonstrate that ε(n) and x2(n) are uncorrelated.

Analyze audio signals

[y,Fs]=audioread('x1.wav');
y=y(:,1);
y_temp=y(1:90000);
dfactor=3;
y=decimate(y_temp,dfactor);
[b,a] = butter(10,70/(Fs/(dfactor*2)),'high');
y=filter(b,a,y);
y=interp(y,dfactor);
sigLength=length(y);
Y = fft(y,sigLength);
Pyy = Y.* conj(Y) / sigLength;
halflength=floor(sigLength/2);
f=Fs*(0:halflength)/sigLength;
figure;plot(f,Pyy(1:halflength+1));xlabel('Frequency(Hz)');
sigLength=length(y_temp);
Y = fft(y_temp,sigLength);
Pyy = Y.* conj(Y) / sigLength;
halflength=floor(sigLength/2);
f=Fs*(0:halflength)/sigLength;
figure;plot(f,Pyy(1:halflength+1));xlabel('Frequency(Hz)');
t=(0:sigLength-1)/Fs;
figure;plot(t,y,t,y_temp);xlabel('Time(s)');
%sound(y,Fs);
sound(y_temp,Fs);

x2.wav is the same.

Remove background noise

clear;
[x,fs]=audioread('x1.wav');
[n,fs]=audioread('x2.wav');
y1=x(1:4096,1);
y2=n(1:4096,1);
Y=fft(y1)-fft(y2);
magY=abs(Y);
b1=[];a=10;b=0.5;
for i=0:126
n=4096;
x1=x(1+n*i:n+n*i);
X1=fft(x1);
magX=abs(X1);
S=(magX.^a-magY.^a);
S1=abs(S).^(1/b);
s1=ifft(S1);
m=mean(s1)*300;
for j=1:4096
if abs(s1(j))>m
s1(j)=s1(j)/4;
end
end
a1=s1';
b1=[b1 a1];
end
x2=b1';
plot(x2);
sound(x2,fs);

Result

Through playback, you can hear clear music.