В этом примере показано, как проектировать фильтры полудиапазона FIR. Полуполосные фильтры широко используются в приложениях многоскоростной обработки сигналов при интерполяции/прореживании в два раза. Полуполосные фильтры эффективно реализуются в многофазной форме, так как приблизительно половина его коэффициентов равна нулю.
Полуполосные фильтры имеют две важные характеристики, пульсации полосы пропускания и полосы останова должны быть одинаковыми, а частоты границы полосы пропускания и границы полосы останова равноудалены от частоты Fs/4 полосы пропускания (Или пи/2 рад/выборка в нормализованной частоте).
firhalfband функция возвращает коэффициенты FIR полуполосного равнополосного фильтра. В качестве простого примера рассмотрим полуполосный фильтр, имеющий отношение к данным, дискретизированным на частоте 96 кГц и частоте полосы пропускания 22 кГц.
Fs = 96e3; Fp = 22e3; N = 100; num = firhalfband(N,Fp/(Fs/2)); fvt = fvtool(num,'Fs',Fs,'Color','white'); fvt.MagnitudeDisplay = 'Zero-phase';
Увеличивая масштаб ответа, можно проверить, что полоса пропускания и полоса останова совпадают. Также существует симметрия относительно точки Fs/4 (24 кГц). Полоса пропускания простирается до 22 кГц, как указано, и полоса останова начинается с 26 кГц. Мы также можем проверить, что каждый другой коэффициент равен нулю, глядя на импульсную характеристику. Это делает фильтр очень эффективным для реализации интерполяции/прореживания с коэффициентом 2.
fvt.Analysis = 'impulse';
dsp.FIRHalfbandInterpolator и dsp.FIRHalfbandDecimator firhalfband функция предоставляет несколько других вариантов конструкции. Однако в большинстве случаев предпочтительно работать непосредственно с dsp.FIRHalfbandInterpolator и dsp.FIRHalfbandDecimator. Эти два системных объекта не только проектируют коэффициенты, но также обеспечивают эффективную реализацию многофазного интерполятора/прореживателя. Они поддерживают фильтрацию данных с плавающей запятой двойной/одинарной точности, а также данных с фиксированной запятой. Они также поддерживают генерацию кода C и HDL, а также оптимизированную генерацию кода ARM ® Cortex ® M и Cortex A.
halfbandInterpolator = dsp.FIRHalfbandInterpolator('SampleRate',Fs, ... 'Specification','Filter order and transition width', ... 'FilterOrder',N,'TransitionWidth',4000); fvtool(halfbandInterpolator,'Fs',2*Fs,'Color','white');
Для выполнения интерполяции dsp.FIRHalfbandInterpolator Используется системный объект. Поскольку это многоскоростной фильтр, важно определить, что подразумевается под частотой дискретизации. Для этого и всех других объектов системы частота дискретизации относится к частоте дискретизации входного сигнала. Однако FVTool определяет частоту дискретизации как скорость, с которой работает фильтр. В случае интерполяции необходимо выполнить повторную выборку, а затем фильтрацию (концептуально), поэтому частота дискретизации FVTool должна быть определена как 2 * Fs из-за увеличения дискретизации на 2.
FrameSize = 256; scope = dsp.SpectrumAnalyzer('SampleRate',2*Fs,'SpectralAverages',5); sine1 = dsp.SineWave('Frequency',10e3','SampleRate',Fs, ... 'SamplesPerFrame',FrameSize); sine2 = dsp.SineWave('Frequency',20e3','SampleRate',Fs, ... 'SamplesPerFrame',FrameSize); tic while toc < 10 x = sine1() + sine2() + 0.01.*randn(FrameSize,1); % 96 kHz y = halfbandInterpolator(x); % 192 kHz scope(y); end release(scope);
Обратите внимание на то, что спектральные реплики ослаблены примерно на 40 дБ, что примерно равно ослаблению, обеспечиваемому фильтром половинной полосы частот. Компенсируя групповую задержку в фильтре, можно построить график наложенных входных и интерполированных выборок. Обратите внимание, что входные выборки остаются неизменными на выходе фильтра. Это происходит потому, что одна из многофазных ветвей полуполосы является чистой ветвью задержки, которая не изменяет входные выборки.
grpDel = 50; n = 0:2:511; stem(n(1:end-grpDel/2),x(1:end-grpDel/2),'k','filled') hold on nu = 0:511; stem(nu(1:end-grpDel),y(grpDel+1:end)) legend('Input samples','Interpolated samples')
В случае прореживания частота выборки, указанная в dsp.FIRHalfbandDecimator соответствует частоте выборки фильтра, так как вы фильтруете, а затем понижаете выборку (концептуально). Так что для прореживателей F, указанные в FVTool, не нужно умножать ни на какой коэффициент.
FrameSize = 256; FsIn = 2*Fs; halfbandDecimator = dsp.FIRHalfbandDecimator('SampleRate',FsIn, ... 'Specification','Filter order and transition width', ... 'FilterOrder',N,'TransitionWidth',4000); fvtool(halfbandDecimator,'Fs',FsIn,'Color','white'); scope = dsp.SpectrumAnalyzer('SampleRate',Fs,'SpectralAverages',5); sine1 = dsp.SineWave('Frequency',10e3','SampleRate',Fs, ... 'SamplesPerFrame',FrameSize); sine2 = dsp.SineWave('Frequency',20e3','SampleRate',Fs, ... 'SamplesPerFrame',FrameSize); tic while toc < 10 x = sine1() + sine2() + 0.01.*randn(FrameSize,1); % 96 kHz y = halfbandInterpolator(x); % 192 kHz xd = halfbandDecimator(y); % 96 kHz scope(xd); end release(scope);
Коэффициенты фильтра могут быть извлечены из интерполятора/прореживателя с помощью tf функция.
num = tf(halfbandInterpolator); % Or num = tf(halfbandDecimator);
Вместо указания порядка фильтра и ширины перехода можно создать фильтр минимального порядка, который обеспечивает заданную ширину перехода, а также заданное затухание полосы останова.
Ast = 80; % 80 dB halfbandInterpolator = dsp.FIRHalfbandInterpolator('SampleRate',Fs, ... 'Specification','Transition width and stopband attenuation', ... 'StopbandAttenuation',Ast,'TransitionWidth',4000); fvtool(halfbandInterpolator,'Fs',2*Fs,'Color','white');
Обратите внимание, что как и для всех интерполяторов, коэффициент усиления полосы пропускания в абсолютных единицах равен коэффициенту интерполяции (2 в случае полубандов). Это соответствует коэффициенту усиления полосы пропускания 6,02 дБ.
Также можно задать порядок фильтрации и затухание полосы останова.
halfbandDecimator = dsp.FIRHalfbandDecimator('SampleRate',Fs, ... 'Specification','Filter order and stopband attenuation', ... 'StopbandAttenuation',Ast,'FilterOrder',N); fvtool(halfbandDecimator,'Fs',Fs,'Color','white');
В отличие от интерполяторов, прореживатели имеют коэффициент усиления 1 (0 дБ) в полосе пропускания.
Полуполосные интерполяторы и прореживатели могут использоваться для эффективной реализации банков фильтров синтеза/анализа. Показанные до сих пор полуполосные фильтры были фильтрами нижних частот. С помощью одного дополнительного сумматора можно получить отклик верхних частот в дополнение к отклику нижних частот и использовать два отклика для реализации банка фильтров.
Следующий код моделирует банк квадратурных зеркальных фильтров (QMF). Сигнал 8 кГц, состоящий из синусоидальных волн 1 кГц и 3 кГц, разделяют на два сигнала 4 кГц, используя полуполосный прореживатель нижних/верхних частот. Сигнал нижних частот сохраняет синусоидальную волну 1 кГц, в то время как сигнал верхних частот сохраняет синусоидальную волну 3 кГц (которая совмещается с 1 кГц после понижающей дискретизации). Затем сигналы объединяются обратно вместе с набором фильтров синтеза с использованием полупериодного интерполятора. Ветвь верхних частот преобразует совмещенную синусоидальную волну 1 кГц обратно в 3 кГц. Интерполированный сигнал имеет частоту дискретизации 8 кГц.
Fs1 = 8000; % Units = Hz Spec = 'Filter order and transition width'; Order = 52; TW = 4.1e2; % Units = Hz % Construct FIR Halfband Interpolator halfbandInterpolator = dsp.FIRHalfbandInterpolator( ... 'Specification',Spec, ... 'FilterOrder',Order, ... 'TransitionWidth',TW, ... 'SampleRate',Fs1/2, ... 'FilterBankInputPort',true); % Construct FIR Halfband Decimator halfbandDecimator = dsp.FIRHalfbandDecimator( ... 'Specification',Spec, ... 'FilterOrder',Order, ... 'TransitionWidth',TW, ... 'SampleRate',Fs1); % Input f1 = 1000; f2 = 3000; InputWave = dsp.SineWave('Frequency',[f1,f2],'SampleRate',Fs1, ... 'SamplesPerFrame',1024,'Amplitude',[1 0.25]); % Construct Spectrum Analyzer object to view the input and output scope = dsp.SpectrumAnalyzer('SampleRate',Fs1, ... 'PlotAsTwoSidedSpectrum',false,'ShowLegend',true,'YLimits', ... [-120 30], ... 'Title', ... 'Input Signal and Output Signal of Quadrature Mirror Filter'); scope.ChannelNames = {'Input','Output'}; tic while toc < 10 Input = sum(InputWave(),2); NoisyInput = Input+(10^-5)*randn(1024,1); [Lowpass,Highpass] = halfbandDecimator(NoisyInput); Output = halfbandInterpolator(Lowpass,Highpass); scope([NoisyInput,Output]); end release(scope);
Все представленные до сих пор конструкции были оптимальными равноудаленными. Используя fdesign.interpolator и fdesign.decimator, доступны другие алгоритмы проектирования.
Fs = 44.1e3; N = 90; TW = 1000/Fs; % Transition width filtSpecs = fdesign.interpolator(2,'halfband','N,TW',N,TW); equirippleHBFilter = design(filtSpecs,'equiripple','SystemObject',true); leastSquaresHBFilter = design(filtSpecs,'firls','SystemObject',true); kaiserHBFilter = design(filtSpecs,'kaiserwin','SystemObject',true);
Можно сравнить проекты с FVTool. Различные конструкции допускают компромиссы между минимальным затуханием полосы останова и более общим затуханием.
fvt = fvtool(equirippleHBFilter,leastSquaresHBFilter,kaiserHBFilter, ... 'Fs',2*Fs,'Color','white'); legend(fvt,'Equiripple design','Least-squares design', ... 'Kaiser-window design')
В качестве альтернативы можно задать порядок и затухание полосы останова. Это обеспечивает компромисс между общим затуханием полосы останова и шириной перехода.
Ast = 60; % Minimum stopband attenuation filtSpecs = fdesign.interpolator(2,'halfband','N,Ast',N,Ast); equirippleHBFilter = design(filtSpecs,'equiripple','SystemObject',true); kaiserHBFilter = design(filtSpecs,'kaiserwin','SystemObject',true); fvt = fvtool(equirippleHBFilter,kaiserHBFilter,'Fs',2*Fs,'Color','white'); legend(fvt,'Equiripple design','Kaiser-window design')
Конструкции окон Kaiser также могут использоваться в дополнение к конструкциям equiripple при проектировании фильтра минимального порядка, необходимого для соответствия спецификациям конструкции. Фактический порядок для конструкции окна Kaiser больше, чем тот, который необходим для конструкции equiripple, но общее затухание полосы останова лучше взамен.
Fs = 44.1e3; TW = 1000/(Fs/2); % Transition width Ast = 60; % 60 dB minimum attenuation in the stopband filtSpecs = fdesign.decimator(2,'halfband','TW,Ast',TW,Ast); equirippleHBFilter = design(filtSpecs,'equiripple','SystemObject',true); kaiserHBFilter = design(filtSpecs,'kaiserwin','SystemObject',true); fvt = fvtool(equirippleHBFilter,kaiserHBFilter,'Fs',Fs,'Color','white'); legend(fvt,'Equiripple Design','Kaiser-window design')
Вместо проектирования оконных фильтров Кайзера можно также получить возрастающее затухание стоп-полосы с помощью модифицированных «эквириптных» конструкций.
equirippleHBFilter1 = design(filtSpecs,'equiripple', ... 'StopbandShape','1/f','StopbandDecay',4,'SystemObject',true); equirippleHBFilter2 = design(filtSpecs,'equiripple', ... 'StopbandShape','linear','StopbandDecay',53.333,'SystemObject',true); fvt = fvtool(equirippleHBFilter1,equirippleHBFilter2, ... 'Fs',Fs,'Color','white'); legend(fvt,'Stopband decaying as (1/f)^4','Stopband decaying linearly')
Фильтр полуполосы верхних частот может быть получен из фильтра полуполосы нижних частот путем изменения знака каждого второго коэффициента. В качестве альтернативы можно непосредственно создать полуполосу верхних частот, установив для свойства «Type» значение «Highpass».
filtSpecs = fdesign.decimator(2,'halfband', ... 'Type','Highpass','TW,Ast',TW,Ast); halfbandHPFilter = design(filtSpecs,'equiripple', ... 'StopbandShape','linear','StopbandDecay',53.333,'SystemObject',true); fvt = fvtool(halfbandHPFilter,equirippleHBFilter2,'Fs',Fs,'Color','white'); legend(fvt,'Highpass halfband filter','Lowpass halfband filter')
