Этот пример показывает, как спроектировать идеальные восстановительные двухканальные фильтрующие банки, также известные как банки квадратурного зеркального фильтра (QMF), поскольку они используют дополнительные фильтры степени.
Часто при обработке цифрового сигнала возникает необходимость разложения сигналов на низкие и высокая частота полосы, после чего необходимо объединить, чтобы восстановить исходный сигнал. Такой пример обнаружен в кодировании поддиапазона (SBC).
Этот пример сначала симулирует идеальный процесс реконструкции путем фильтрации сигнала, составленного из дельт Кронекера. Предусмотрены графики входного, выходного и сигнала ошибки, а также спектр величин передаточной функции всей системы. Эффективность идеальной реконструкции показана через этот банк фильтров. После этого в примере приложения показано, как два поддиапазона аудио файла могут обрабатываться по-разному без особых эффектов на реконструкцию.
Идеальная реконструкция - это процесс, посредством которого сигнал полностью восстанавливается после разделения на его низкие частоты и высокие частоты. Ниже приведен блок идеального процесса реконструкции, который использует идеальные фильтры. Идеальный процесс реконструкции требует четырех фильтров, два lowpass (H0 и G0) и двух highpass-фильтров (H1 и G1). В сложение, это требует понижающей дискретизации и повышающей дискретизации между два lowpass и между двумя высокочастотными фильтрами. Обратите внимание, что мы должны принять во внимание тот факт, что наши выходные фильтры должны иметь коэффициент усиления в два для компенсации предыдущего усилителя.
DSP System Toolbox™ обеспечивает специализированную функцию, называемую FIRPR2CHFB, для разработки четырех фильтров, необходимых для реализации идеальной реконструкции конечной импульсной характеристики двухканального банка фильтров, как описано выше. FIRPR2CHFB проектирует четыре конечные импульсные характеристики для секций анализа (H0 и H1) и синтеза (G0 и G1) двухканального идеального набора восстановительных фильтров. Это проект соответствует так называемым ортогональным блокам фильтров, также известным как силосимметричные блоки фильтров, которые требуются в порядок для достижения совершенной реконструкции.
Спроектируем банк фильтров с фильтрами порядка 99 и ребра полосы пропускания lowpass и highpass фильтров 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');
Обратите внимание, что путь анализа состоит из фильтра, за которым следует понижающий усилитель, который является дециматором, и путь синтеза состоит из повышающего усилителя, за которым следует фильтр, который является интерполятором. Для реализации DSP System Toolbox™ Системных объектов предоставляет две 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
Во-первых, создайте Системный объект для загрузки и игры аудио файла.
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
Затем мы хотим сбросить фильтры поддиапазона из вышеописанного примера идеальной реконструкции, чтобы мы могли повторно использовать их. Метод Release вызывается на графиках, чтобы иметь возможность изменять определенные свойства.
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. Можно видеть, что коэффициент усиления отличается для низких частот, чем для более высоких таковых. Однако, услышав воспроизведение восстановленного аудиосигнала, вы бы заметили, что человеческое ухо не слишком восприимчиво к изменению разрешения, более того, в случае высоких частот, где длина слова была ещё меньше.