Этот пример демонстрирует БПФ масштаба, который является методом обработки сигналов, используемым для анализа фрагмента спектра с высоким разрешением. DSP System Toolbox™ предлагает эту функциональность в MATLAB™ через dsp.ZoomFFT Системный объект и в Simulink через библиотечный блок FFT с масштабированием.
Разрешение сигнала ограничено его длиной. Мы проиллюстрируем этот факт простым примером. Рассмотрим сигнал, формируемый суммой двух синусоид:
L = 32; % Frame size Fs = 128; % Sample rate res = Fs/L; % Frequency resolution f1 = 40; % First sine wave frequency f2 = f1 + res; % Second sine wave frequency sn1 = dsp.SineWave('Frequency',f1,'SampleRate',Fs,'SamplesPerFrame',L); sn2 = dsp.SineWave('Frequency',f2,'SampleRate',Fs,'SamplesPerFrame',L,... 'Amplitude',2); x = sn1() + sn2();
Давайте вычислим БПФ x и постройте величину БПФ. Обратите внимание, что две синусоиды находятся в смежных выборках БПФ. Это означает, что вы не можете различать частоты ближе, чем Fs/L.
X = fft(x); stem(Fs/L*(0:length(X)-1)-Fs/2,abs(fftshift(X))/L) grid on; xlabel('Frequency (Hz)') ylabel('Magnitude') title('Two-sided spectrum')
Предположим, у вас есть приложение, для которого вас интересует только поддиапазон интервала Найквиста. Идея zoom FFT состоит в том, чтобы сохранить то же разрешение, которого вы достигли бы с полноразмерным БПФ на вашем исходном сигнале, вычисляя небольшой БПФ на более коротком сигнале. Более короткий сигнал приходит от децимирования исходного сигнала. Экономия за счет вычисления гораздо более короткого БПФ при достижении того же разрешения. Это интуитивно понятно: для коэффициента десятикратного уменьшения D новая частота дискретизации - Fsd = Fs/D, а новый формат кадра (и длина БПФ) - Ld = L/D, поэтому разрешение децимированного сигнала - Fsd/Ld = Fs/L.
DSP System Toolbox предлагает функциональность БПФ для MATLAB и Simulink, через dsp. Системный объект и библиотечный блок БПФ, соответственно. В следующих разделах будет обсуждаться алгоритм БПФ более подробно.
Перед обсуждением алгоритма, используемого в dsp. БПФ, мы представляем миксер подход, который является популярным методом БПФ.
Для примера здесь предположим, что вас интересует только интервал [16 Гц, 48 Гц].
BWOfInterest = 48 - 16;
Fc = (16 + 48)/2; % Center frequency of bandwidth of interest
Достижимым фактором десятикратного уменьшения является:
BWFactor = floor(Fs/BWOfInterest);
Подход миксера состоит сначала в смещении интересующей полосы вниз к DC с помощью миксера, а затем в выполнении lowpass фильтрации и децимации в множителе BWFactor (с использованием эффективной полифазной КИХ-децимирующей структуры).
Давайте спроектируем коэффициенты децимирующего фильтра с помощью функции designMultirateFIR:
B = designMultirateFIR(1,BWFactor); D = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',B);
Теперь давайте смешаем сигнал до DC и фильтруем его через дециматор конечную импульсную характеристику:
for k = 1:10 % Run a few times to eliminate transient in filter x = sn1()+sn2(); % Mix down to DC indVect = (0:numel(x)-1).' + (k-1) * size(x,1); y = x .* exp(-2*pi*indVect*Fc*1j/Fs); % Filter through FIR decimator xd = D(y); end
Теперь возьмем БПФ фильтрованного сигнала (обратите внимание, что длина БПФ уменьшается на BWFactor, или на длину десятикратного уменьшения, по сравнению с обычной БПФ, при сохранении того же разрешения):
fftlen = numel(xd); Xd = fft(xd); figure Ld = L/BWFactor; Fsd = Fs/BWFactor; F = Fc + Fsd/fftlen*(0:fftlen-1)-Fsd/2; stem(F,abs(fftshift(Xd))/Ld) grid on xlabel('Frequency (Hz)') ylabel('Magnitude') title('Zoom FFT Spectrum. Mixer Approach.')
Комплексный смеситель добавляет дополнительное умножение для каждой высокоскоростной выборки, которая не является эффективной. в следующем разделе представлен альтернативный, более эффективный, масштабируемый подход БПФ.
Альтернативный метод БПФ использует в своих интересах известный результат полосно-пропускающей фильтрации (также иногда называемой недостаточной дискретизацией): Предположим, что мы заинтересованы в полосе [F1,F2] сигнала с частотой дискретизации Fs Гц. Если мы передадим сигнал через комплексный (односторонний) полосовой фильтр с центром в Fc = (F1 + F2 )/2 и с шириной полосы BW = F2 - F1, а затем понизим его в множитель D = floor (Fs/BW), мы сбросим нужную полосу до основной полосы.
В целом, если Fc не может быть выражен в форме k * Fs/D (где K - целое число), то сдвинутый, децимированный спектр не будет центрироваться при DC. Фактически центральная частота Fc будет переведена в [2]:
Fd = Fc - (Fs/D) * floor ((D * Fc + Fs/2 )/Fs)
В этом случае мы можем использовать миксер (вращающийся с низкой частотой дискретизации децимированного сигнала), чтобы центрировать нужную полосу до нуля Hertz.
Используя пример из предыдущего раздела, получаем коэффициенты комплексного полосно-пропускающего фильтра путем модулирования коэффициентов проектируемого lowpass:
Bbp = B .*exp(1j*(Fc / Fs)*2*pi*(0:numel(B)-1)); D = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',Bbp);
Теперь выполните фильтрацию и БПФ:
for k = 1:10 % Run a few times to eliminate transient in filter x = sn1()+sn2(); xd = D(x); end Xd = fft(xd); figure stem(F,abs(fftshift(Xd))/Ld) grid on xlabel('Frequency (Hz)') ylabel('Magnitude') title('Zoom FFT Spectrum. Bandpass Sampling Approach.')
Фильтр из предыдущего раздела является однокаскадным. Мы можем достичь менее вычислительно сложного фильтра, используя вместо этого многоступенчатый проект. Это подход, используемый в dsp. ZoomFFT.
Рассмотрим следующий пример, где входная частота выборки составляет 48 кГц, а интересующая полоса пропускания является интервалом [1500,2500] Гц. Достижимый коэффициент десятикратного уменьшения равен этажу (48000/1000) = 48.
Давайте сначала спроектируем одноступенчатый дециматор:
Fs = 48e3; Fc = 2000; % Bandpass filter center frequency BW = 1e3; % Bandwidth of interest Ast = 80; % Stopband attenuation BWFactor = floor(Fs/BW); B = designMultirateFIR(1,BWFactor,12,80); Bbp = B .*exp(1j*(Fc / Fs)*2*pi*(0:numel(B)-1)); D_single_stage = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',Bbp);
Теперь давайте спроектируем фильтр с помощью многоступенчатого проекта, сохраняя при этом то же затухание в полосе задерживания и полосу пропускания перехода, что и одноступенчатый случай (см kaiserord для получения дополнительной информации об расчете ширины перехода):
tw = (Ast - 7.95) / ( numel(B) * 2.285); fmd = fdesign.decimator(BWFactor,'Nyquist',BWFactor,... 'TW,Ast', tw * Fs / (2*pi),Ast,Fs); D_multi_stage = multistage(fmd,'HalfbandDesignMethod','equiripple','SystemObject',true); fprintf('Number of filter stages: %d\n', getNumStages(D_multi_stage) ); fprintf('Stage 1: Decimation factor = %d , filter length = %d\n',... D_multi_stage.Stage1.DecimationFactor,... numel(D_multi_stage.Stage1.Numerator)); fprintf('Stage 2: Decimation factor = %d , filter length = %d\n',... D_multi_stage.Stage2.DecimationFactor,... numel(D_multi_stage.Stage2.Numerator)); fprintf('Stage 3: Decimation factor = %d , filter length = %d\n',... D_multi_stage.Stage3.DecimationFactor,... numel(D_multi_stage.Stage3.Numerator));
Number of filter stages: 3 Stage 1: Decimation factor = 6 , filter length = 33 Stage 2: Decimation factor = 2 , filter length = 11 Stage 3: Decimation factor = 4 , filter length = 101
Обратите внимание, что D_multi_stage является трехступенчатым многоразовым lowpass. Мы преобразуем его в полосно-пропускающий фильтр, выполняя сдвиг частоты на коэффициентах каждого каскада, принимая при этом в расчет совокупный коэффициент десятикратного уменьшения:
num1 = D_multi_stage.Stage1.Numerator; num2 = D_multi_stage.Stage2.Numerator; num3 = D_multi_stage.Stage3.Numerator; N1 = length(num1)-1; N2 = length(num2)-1; N3 = length(num3)-1; % Decimation factor between stage 1 & 2: M12 = D_multi_stage.Stage1.DecimationFactor; % Decimation factor between stage 2 & 3: M23 = D_multi_stage.Stage2.DecimationFactor; num1 = num1 .*exp(1j*(Fc / Fs)*2*pi*(0:N1)); num2 = num2 .*exp(1j*(Fc * M12 / Fs)*2*pi*(0:N2)); num3 = num3 .*exp(1j*(Fc * M12 * M23 / Fs)*2*pi*(0:N3)); D_multi_stage.Stage1.Numerator = num1; D_multi_stage.Stage2.Numerator = num2; D_multi_stage.Stage3.Numerator = num3;
Сравним стоимость одноступенчатого и многоступенчатого фильтров:
fprintf('Single-stage filter cost:\n\n') cost(D_single_stage) fprintf('Multistage filter cost:\n') cost(D_multi_stage)
Single-stage filter cost:
ans =
struct with fields:
NumCoefficients: 1129
NumStates: 1104
MultiplicationsPerInputSample: 23.5208
AdditionsPerInputSample: 23.5000
Multistage filter cost:
ans =
struct with fields:
NumCoefficients: 113
NumStates: 140
MultiplicationsPerInputSample: 7.0208
AdditionsPerInputSample: 6.7500
Давайте также сравним частотную характеристику двух фильтров:
vis = fvtool(D_single_stage,D_multi_stage,'DesignMask','off','legend','on'); legend(vis,'Single-stage','Multistage')
Давайте будем использовать многоступенчатый фильтр для выполнения масштабирования БПФ:
fftlen = 32; L = BWFactor * fftlen; tones = [1625 2000 2125]; % sine wave tones sn = dsp.SineWave('SampleRate',Fs,'Frequency',tones,... 'SamplesPerFrame',L); Fsd = Fs / BWFactor; % Frequency points at which FFT is computed F = Fc + Fsd/fftlen*(0:fftlen-1)-Fsd/2; ap = dsp.ArrayPlot('XDataMode','Custom',... 'CustomXData',F,... 'YLabel','Magnitude',... 'XLabel','Frequency (Hz)',... 'YLimits',[0 1],... 'Title',sprintf('Zoom FFT. Resolution = %f Hz',(Fs/BWFactor)/fftlen)); for k=1:100 x = sum(sn(),2) + 1e-2 * randn(L,1); y = D_multi_stage(x); z = fft(y,fftlen,1); z = fftshift(z); ap( abs(z)/fftlen ) end release(ap)
dsp.ZoomFFT является системным объектом, который реализует масштабирование БПФ на основе многоступенчатого многоступенчатого полосно-пропускающего фильтра, выделенного в предыдущем разделе. Вы задаете необходимую центральную частоту и коэффициент десятикратного уменьшения, и dsp. ZoomFFT разработает фильтр и применит его к входному сигналу.
Давайте использовать dsp. ZoomFFT для масштабирования синусоидальных тонов из примера предыдущего раздела:
zfft = dsp.ZoomFFT(BWFactor,Fc,Fs); for k=1:100 x = sum(sn(),2) + 1e-2 * randn(L,1); z = zfft(x); z = fftshift(z); ap( abs(z)/fftlen) end release(ap)
Блок БПФ приносит функциональность dsp. Масштабирование FFT до Simulink. В модели dspzoomfft мы используем блок БПФ для просмотра полосы частот [800 Гц, 1600 Гц] входного сигнала, дискретизированного с частотой 44100 Гц.
[1] Многоуровневая обработка сигналов - Харрис (Prentice Hall).
[2] Вычисление переведенных частот в оцифровке и понижающей дискретизации аналоговой полосы пропускания - Lyons (https://www.dsprelated.com/showarticle/523.php)