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

Вход является звуковым сигналом, произведенным на уровне 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. Установите длину БПФ на 1 024.

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 одна длина транзитного участка за один раз. Примените предполагаемое усиление к преобразованному сигналу. Восстановите denoised версию исходного речевого сигнала путем выполнения обратного преобразования Фурье на отдельных диапазонах частот. Проигрывайте denoised звуковой сигнал к аудио устройству компьютера.

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 если комбинация входных параметров совместима КОЛОЙ и логический false если нет. method аргумент функции установлен в 'ola' или '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 со взвешенным перекрытием - добавляют (ВОЛЯ)

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

iscola Функция

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

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

tfWOLA = logical
   1

Восстановите данные с ВОЛЕЙ

Выпустите 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 из различия между входным сигналом и восстановленным сигналом очень мал, указывая, что сигнал был восстановлен отлично.