ee2800/rlsdsp.m at master · ARandomDev99/ee2800 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
% selection of noise suppression mode
p = 1;

Fs = 44100;             % Samplinng frequrncy
f1 = 1000;              % Tonal frequencies present in the noise
w = 2 * pi * f1 / Fs;   % Normalized tonal frequencies

external_noise = load('Ameen Aazam - external_noise.txt');
noisy_speech = load('Ameen Aazam - noisy_speech.txt');
clean_speech = load('Ameen Aazam - clean_speech.txt');

f = (0:length(external_noise)-1) * Fs / length(external_noise);     % Frequency axis (for viewing fft if needed)
x = external_noise;
y = zeros(length(noisy_speech), length(f1) + 1);
y(:, 1) =  x;
r = 0.999;
s = clean_speech;

N = 2;
lambda = 1;
delta = 0.1;
z = noisy_speech;

wrls = zeros(N, 1);
P = eye(N) / delta;
erls = zeros(length(z), 1);
yrls = zeros(length(z), 1);

for i = N:length(noisy_speech)
    if (p == 1)
        for j = 1:length(f1)
            %fprintf("%d\n", i);
            if (i == 1)
                y(i, j + 1) = y(i, j);
            elseif (i == 2)
                y(i, j + 1) = y(i, j) - 2 * cos(w(j)) * y(i - 1, j) + 2 * r * cos(w(j)) * y(i - 1, j + 1);
            else
                y(i, j + 1) = y(i, j) - 2 * cos(w(j)) * y(i - 1, j) + y(i - 2, j) + 2 * r * cos(w(j)) * y(i - 1, j + 1) - r^2 * y(i - 2, j + 1);
            end
        end
        yi = y(i:-1:i-N+1, end);
        yrls(i) = wrls' * yi;
        erls(i) = z(i) - yrls(i);                   % Output for partial suppression
        k = (P * yi) / (lambda + yi' * P * yi);
        wrls = wrls + k * erls(i);
        P = (P - k * yi' * P) / lambda;
    else
        %fprintf("%d\n", i);
        yi = y(i:-1:i-N+1, 1);
        yrls(i) = wrls' * yi;
        erls(i) = z(i) - yrls(i);                   % Output for full suppression
        k = (P * yi) / (lambda + yi' * P * yi);
        wrls = wrls + k * erls(i);
        P = (P - k * yi' * P) / lambda;
    end
end

% Calculation of SNR for full suppression
if (p == 0)
    ps = 0;
    pn = 0;
    for i = 1 : length(noisy_speech) - N + 1
        ps = ps + s(i) ^ 2;
        pn = pn + (erls(i) - s(i)) ^ 2;
    end

    SNR = 10 * log10(ps/pn);

% Calculation of SNR for partial suppression (Separate SNR for the signal and the tone(s))
else
    t = 0:1/Fs:(length(x) - 1)/Fs;
    l = zeros(length(x), 1);
    % Building the input signal
    for i = 1:length(f1)
        l = l + (sin(2 * pi * f1(i) * t))';
    end

    % Separating the tonal component(s) from the output (Leaving the signal an the random noise)
    Nsnr = 5;
    lambdasnr = 1;
    deltasnr = 0.1;
    zsnr = erls;

    wsnr = zeros(Nsnr, 1);
    Psnr = eye(Nsnr) / deltasnr;
    esnr = zeros(length(l), 1);
    ysnr = zeros(length(l), 1);
    for i = Nsnr:length(l)
        li = l(i:-1:i-Nsnr+1);
        ysnr(i) = wsnr' * li;
        esnr(i) = zsnr(i) - ysnr(i);
        ksnr = (Psnr * li) / (lambdasnr + li' * Psnr * li);
        wsnr = wsnr + ksnr * esnr(i);
        Psnr = (Psnr - ksnr * li' * Psnr) / lambdasnr;
    end
    ps = 0;
    pt = 0;
    pn = 0;
    for i = 1 : length(noisy_speech) - N + 1
        ps = ps + s(i) ^ 2;
        pt = pt + ysnr(i) ^ 2;
        pn = pn + (esnr(i) - s(i)) ^ 2;
    end

    SNRS = 10 * log10(ps/pn);           % SNR of the signal
    SNRT = 10 * log10(pt/pn);           % SNR of the tone(s)
end