Шумный входной сигнал

Вход представляет собой звуковой сигнал, дискретизированный на частоте 22 050 Гц. dsp.AudioFileReader объект считывает этот сигнал в кадрах 512 выборок. Звуковой сигнал искажается белым гауссовым шумом, который имеет стандартное отклонение 0,05. Используйте audioDeviceWriter объект для воспроизведения шумного звукового сигнала на звуковое устройство компьютера.

FrameLength = 512;
afr = dsp.AudioFileReader('speech_dft.wav',...
    'SamplesPerFrame',FrameLength);
adw = audioDeviceWriter('SampleRate',afr.SampleRate);

noiseStd = 0.05;
while ~isDone(afr)
    cleanAudio = afr();
    noisyAudio = cleanAudio + noiseStd * randn(FrameLength,1);
    adw(noisyAudio);
end
reset(afr)

Инициализация объектов преобразования Фурье с коротким и обратным коротким временем

Инициализируйте dsp.STFT и dsp.ISTFT объекты. Установите длину окна равной длине входного кадра, а длину транзитного участка - 16. Длина перекрытия представляет собой разницу между длиной окна и длиной перехода, OL = WL-HL. Установите длину БПФ 1024.

WindowLength = FrameLength;
HopLength = 16;
numHopsPerFrame = FrameLength / 16;
FFTLength = 1024;

Окно, используемое для вычисления STFT и ISTFT, является периодическим hamming окно длиной 512. ConjugateSymmetricInput флаг istf объект имеет значение true, указывая, что выходные данные istf объект является сопряженно-симметричным сигналом.

win = hamming(WindowLength,'periodic');
stf = dsp.STFT(win,WindowLength-HopLength,FFTLength);
istf = dsp.ISTFT(win,WindowLength-HopLength,1,0);

Оценщик усиления

Следующим шагом является определение параметров оценщика усиления. Этот коэффициент усиления применяется к шумовому спектру для ослабления поддиапазонов с более высокой мощностью шума и подъема поддиапазонов с меньшей мощностью шума.

dec = 16;
alpha = 15;
stftNorm = (sum(win.*win) / dec).^2;

Спектральное ослабление

Подать звуковой сигнал на stf по одному транзитному участку за один раз. Примените расчетный коэффициент усиления к преобразованному сигналу. Реконструируют деноизированную версию исходного речевого сигнала путем выполнения обратного преобразования Фурье в отдельных полосах частот. Воспроизведение запрещенного звукового сигнала на звуковое устройство компьютера.

while ~isDone(afr)
    cleanAudio =  afr();
    noisyAudio = cleanAudio + noiseStd * randn(FrameLength,1);
    y = zeros(FrameLength,1); % y holds the denoised audio frame
    
    % Feed audio to stft one hop-length at a time
    for index = 1:numHopsPerFrame        
        X = stf(noisyAudio((index-1)*HopLength+1:index*HopLength));        
        % Gain estimator
        Z = abs(X).^2 / (noiseStd^2 * alpha) / stftNorm;
        Z(Z<=1) = 1;
        Z = 1 - 1./Z;
        Z = sign(Z) .* sqrt(abs(Z));
        X = X .* Z;        
        % Convert back to time-domain
        y((index-1)*HopLength+1:index*HopLength) = istf(X);        
    end    
    % Listen to denoised audio:
    adw(y);
end

iscola Функция

iscola функция проверяет, чтобы указанное окно и перекрытие удовлетворяли ограничению COLA, чтобы гарантировать, что обратное кратковременное преобразование Фурье (ISTFT) приводит к идеальной реконструкции для немодифицированных спектров. Функция возвращает логическое значение true если комбинация входных параметров совместима с COLA и является логической false если нет. method аргумент функции имеет значение 'ola' или 'wola' в зависимости от того, использует ли метод инверсии взвешенное наложение-сложение (WOLA).

Проверить, hann() окно длиной 120 выборок и длиной перекрытия 60 выборок соответствует COLA.

winLen = 120;
overlapLen = 60;
win = hann(winLen,'periodic');
tf = iscola(win,overlapLen,'ola')

tf = logical
   1

Инициализация

Инициализируйте dsp.STFT и dsp.ISTFT Системные объекты с этим hann , которое совместимо с COLA. Установите длину БПФ равной длине окна.

frameLen = winLen-overlapLen;
stf = dsp.STFT('Window',win,'OverlapLength',overlapLen,'FFTLength',winLen);
istf = dsp.ISTFT('Window',win,'OverlapLength',overlapLen,'WeightedOverlapAdd',0);

Реконструировать данные

Вычислите STFT случайного сигнала. Установите длину входного сигнала равной длине перехода (длина окна - длина перекрытия). Поскольку окно совместимо с COLA, ISTFT этого немодифицированного спектра прекрасно восстанавливает исходный сигнал временной области.

Для подтверждения сравните входной сигнал x с восстановленным выходным сигналом y. Благодаря задержке, введенной объектами, восстановленный выходной сигнал сдвигается во времени по сравнению с входным сигналом. Поэтому для сравнения возьмем норму разности между восстановленным выходом y и предыдущим входом xprev. Норма очень мала, указывая, что выходной сигнал является идеально реконструированным вариантом входного сигнала.

n = zeros(1,100);
xprev = 0;
for i = 1:100
    x = randn(frameLen,1);
    X = stf(x);
    y = istf(X);
    n(1,i) = norm(y-xprev);
    xprev = x;
end       
max(abs(n))

ans = 1.6972e-13

ISTFT с взвешенным перекрытием - добавление (WOLA)

В WOLA второе окно, называемое окном синтеза, $\mathit{f}\mathit{(}n$), применяется после операции IFFT и до наложения-добавления. Окна синтеза и анализа обычно идентичны и обычно получаются путем взятия квадратного корня окон, удовлетворяющих COLA (тем самым обеспечивая идеальную реконструкцию ).

iscola Функция

Проверить, sqrt(hann()) окно длиной 120 выборок и длиной перекрытия 60 выборок соответствует WOLA. Установите method аргумент iscola функция для 'wola'. Выходные данные iscola функция 1 указывает, что это окно совместимо с WOLA.

winWOLA = sqrt(hann(winLen,'periodic'));
tfWOLA = iscola(winWOLA,overlapLen,'wola')

tfWOLA = logical
   1

Реконструкция данных с помощью WOLA

Отпустите клавишу dsp.STFT и dsp.ISTFT Системные объекты и установите для окна значение sqrt(hann(winLen,'periodic')) окно. Использование взвешенного перекрытия - добавление в ISTFT сторона, установите 'WeightedOverlapAdd' кому true.

release(stf);
release(istf);
stf.Window = winWOLA;
istf.Window = winWOLA;
istf.WeightedOverlapAdd = true;

n = zeros(1,100);
xprev = 0;
for i = 1:100
    x = randn(frameLen,1);
    X = stf(x);
    y = istf(X);
    n(1,i) = norm(y-xprev);
    xprev = x;
end       
max(abs(n))

ans = 4.6664e-15

norm разность между входным сигналом и восстановленным сигналом очень мала, что указывает на то, что сигнал восстановлен идеально.