Этот пример показывает, как спроектировать идеальную реконструкцию двухканальных банков фильтров, также известных как банки квадратурного зеркального фильтра (QMF), поскольку они используют дополняющие по мощности фильтры.
Часто при цифровой обработке сигнала возникает необходимость разложения сигналов на низкочастотные и высокочастотные диапазоны, после чего необходимо объединить для восстановления исходного сигнала. Такой пример находится в поддиапазонном кодировании (SBC).
Этот пример сначала смоделирует идеальный процесс реконструкции путем фильтрации сигнала, составленного из дельт Кронекера. Предусмотрены графики входного, выходного и ошибочного сигнала, а также амплитудный спектр передаточной функции всей системы. Эффективность идеальной реконструкции показана через этот банк фильтров. После этого примерное приложение показывает, как два поддиапазона аудиофайла могут обрабатываться по-разному без особого влияния на реконструкцию.
Идеальная реконструкция - это процесс, с помощью которого сигнал полностью восстанавливается после разделения на его низкие частоты и высокие частоты. Ниже приведена блок-схема идеального процесса реконструкции, в котором используются идеальные фильтры. Идеальный процесс реконструкции требует четырех фильтров, двух фильтров нижних частот (H0 и G0) и двух фильтров верхних частот (H1 и G1). Кроме того, требуется устройство понижающей дискретизации и устройство повышающей дискретизации между двумя фильтрами нижних частот и между двумя фильтрами верхних частот. Заметим, что мы должны учитывать тот факт, что наши выходные фильтры должны иметь коэффициент усиления два, чтобы компенсировать предшествующий повышающий дискретизатор.
Система DSP Toolbox™ предоставляет специализированную функцию, называемую FIRPR2CHFB, для разработки четырех фильтров, необходимых для реализации набора двухканальных фильтров с идеальной реконструкцией КИХ, как описано выше. FIRPR2CHFB проектирует четыре фильтра FIR для секций анализа (H0 и H1) и синтеза (G0 и G1) двухканального идеального набора фильтров реконструкции. Конструкция соответствует так называемым ортогональным блокам фильтров, также известным как силосимметричные блоки фильтров, которые требуются для достижения совершенной реконструкции.
Создадим набор фильтров с фильтрами порядка 99 и краями полосы пропускания фильтров нижних и верхних частот 0,45 и 0,55 соответственно:
N = 99; [LPAnalysis, HPAnalysis, LPSynthsis, HPSynthesis] = firpr2chfb(N, 0.45);
Амплитудный отклик этих фильтров показан ниже:
fvt = fvtool(LPAnalysis,1, HPAnalysis,1, LPSynthsis,1, HPSynthesis,1); fvt.Color = [1,1,1]; legend(fvt,'Hlp Lowpass Decimator','Hhp Highpass Decimator',... 'Glp Lowpass Interpolator','Ghp Highpass Interpolator');
Заметим, что путь анализа состоит из фильтра, за которым следует понижающая дискретизация, которая является прореживающим устройством, а путь синтеза состоит из повышающей дискретизации, за которой следует фильтр, который является интерполятором. Системный Toolbox™ DSP предоставляет два системных объекта для реализации этого - dsp.SubbandAnalysisFilter для анализа и dsp.SubbandSynthesisFilter для раздела синтеза.
analysisFilter = dsp.SubbandAnalysisFilter(LPAnalysis, HPAnalysis);
% Analysis section
synthFilter = dsp.SubbandSynthesisFilter(LPSynthsis, HPSynthesis);
% Synthesis section
В качестве примера, пусть p [n] обозначает сигнал
и пусть сигнал x [n] определяется
ПРИМЕЧАНИЕ: Так как MATLAB ® использует индексирование на основе одной единицы, дельта [n] = 1, когда n = 1.
x = zeros(50,1); x(1:3) = 1; x(8:10) = 2; x(16:18) = 3; x(24:26) = 4; x(32:34) = 3; x(40:42) = 2; x(48:50) = 1; sigsource = dsp.SignalSource('SignalEndAction', 'Cyclic repetition',... 'SamplesPerFrame', 50); sigsource.Signal = x;
Для просмотра результатов моделирования потребуется три области - первая для сравнения входного сигнала с восстановленным выходом, вторая для измерения погрешности между двумя и третьей для построения графика амплитудной характеристики всей системы.
% Scope to compare input signal with reconstructed output sigcompare = dsp.ArrayPlot('NumInputPorts', 2, 'ShowLegend', true,... 'Title', 'Input (channel 1) and reconstructed (channel 2) signals'); % Scope to plot the RMS error between the input and reconstructed signals errorPlot = timescope('Title', 'RMS Error', 'SampleRate', 1, ... 'TimeUnits', 'Seconds', 'YLimits', [-0.5 2],... 'TimeSpanSource', 'property','TimeSpan',100,... 'TimeSpanOverrunAction','Scroll'); % To calculate the transfer function of the cascade of Analysis and % Synthesis subband filters tfestimate = dsp.TransferFunctionEstimator('FrequencyRange','centered',... 'SpectralAverages', 50); % Scope to plot the magnitude response of the estimated transfer function tfplot = dsp.ArrayPlot('PlotType','Line', ... 'YLabel', 'Frequency Response (dB)',... 'Title','Transfer function of complete system',... 'XOffset',-25, 'XLabel','Frequency (Hz)');
Теперь мы пропускаем входной сигнал через поддиапазонные фильтры и восстанавливаем выходной сигнал. Результаты выводятся на печать в созданных областях.
for i = 1:100 input = sigsource(); [hi, lo] = analysisFilter(input); % Analysis reconstructed = synthFilter(hi, lo); % Synthesis % Compare signals. Delay input so that it aligns with the filtered % output. Delay is due to the filters. sigcompare(input(2:end), reconstructed(1:end-1)); % Plot error between signals err = rms(input(2:end) - reconstructed(1:end-1)); errorPlot(err); % Estimate transfer function of cascade Txy = tfestimate(input(2:end), reconstructed(1:end-1)); tfplot(20*log10(abs(Txy))); end release(errorPlot);
Мы можем видеть на первых двух графиках, что наш идеальный банк двухканальных фильтров реконструкции полностью реконструировал наш исходный сигнал x [n]. Начальная ошибка вызвана задержкой в фильтрах. Третий график показывает, что каскад фильтров поддиапазонов не изменяет частотные характеристики сигнала.
Поддиапазонные фильтры позволяют обрабатывать высокие частоты в сигнале способом, отличным от способа обработки низких частот. В качестве примера, мы загрузим аудиофайл и квантуем его низкие частоты с длиной слова, которая выше, чем у высоких частот. Реконструкция тогда может быть не идеальной, как выше, но она все равно будет достаточно точной.
Примечание.Для этого требуется лицензия для Fixed-Point Designer
Сначала создайте объект System для загрузки и воспроизведения аудиофайла.
audioInput = dsp.AudioFileReader;
audioWriter = audioDeviceWriter('SampleRate',audioInput.SampleRate);
Чтобы иметь меру отсчета, мы можем воспроизвести исходный звук один раз
while ~isDone(audioInput) input = audioInput(); % Load a frame audioWriter(input); % Play the frame end % Wait until audio is played to the end pause(10*audioInput.SamplesPerFrame/audioInput.SampleRate); reset(audioInput); % Reset to beginning of the file release(audioInput); % Close input file release(audioWriter); % Close audio output device
Далее, мы хотим сбросить фильтры поддиапазонов из приведенного выше примера идеальной реконструкции, чтобы мы могли использовать их повторно. Метод деблокирования вызывается на графиках для изменения определенных свойств.
hide(errorPlot) hide(tfplot) reset(analysisFilter); release(analysisFilter); % Release to change input data size reset(synthFilter); release(synthFilter); % Release to change input data size % Plot for error errorPlot.SampleRate = audioInput.SampleRate/audioInput.SamplesPerFrame; errorPlot.TimeSpan = 5.5; clear sigcompare % Do not need a plot to compare signals reset(tfestimate); % Transfer function estimate release(tfestimate); release(tfplot); % Plot for transfer function estimate tfplot.YLimits = [-20, 60]; tfplot.XOffset = -audioInput.SampleRate/2; tfplot.SampleIncrement = ... audioInput.SampleRate/audioInput.SamplesPerFrame;
Цикл моделирования очень похож на цикл для идеального примера реконструкции в начале. Ниже приведены следующие изменения:
Квантование выполняется для низкочастотного компонента, имеющего 8 битов, и высокочастотного компонента, имеющего 4 бита длины слова.
Восстановленный звук воспроизводится, чтобы позволить пользователю услышать его и заметить любые изменения из входного аудиофайла. Следует отметить, что квантованные поддиапазоны сохраняются в double контейнер, потому что dsp.SubbandSynthesisFilter объект должен иметь один и тот же числовой тип.
show(tfplot) show(errorPlot) while ~isDone(audioInput) input = audioInput(); % Load a frame of audio [hi, lo] = analysisFilter(input); % Analysis QuantizedHi = double(fi(hi, 1, 4, 9)); % Quantize to 4 bits QuantizedLo = double(fi(lo, 1, 8, 8)); % Quantize to 8 bits reconstructed = synthFilter(QuantizedHi, QuantizedLo); % Synthesis % Plot error between signals err = rms(input(2:end) - reconstructed(1:end-1)); errorPlot(err); % Play the reconstructed audio frame audioWriter(reconstructed); % Estimate transfer function of cascade Txy = tfestimate(input(2:end), reconstructed(1:end-1)); tfplot(20*log10(abs(Txy))); end release(errorPlot); % Wait until audio is played to the end pause(10*audioInput.SamplesPerFrame/audioInput.SampleRate); release(audioWriter); % Close audio output device reset(audioInput); % Reset to beginning of the file release(audioInput); % Close input file
Как показывает график ошибок, реконструкция не является совершенной из-за квантования. Также, в отличие от предыдущего случая, оценка передаточной функции полной системы также не 0dB. Коэффициент усиления для низких частот отличается от коэффициента усиления для более высоких частот. Однако, услышав воспроизведение восстановленного аудиосигнала, вы бы заметили, что человеческое ухо не слишком заметно для изменения разрешения, более того, в случае высоких частот, где длина слова была еще меньше.