Этот пример демонстрирует БПФ изменения масштаба, который является методом обработки сигналов, используемым, чтобы анализировать фрагмент спектра в высоком разрешении. 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')
Предположим, что у вас есть приложение, для которого вы только интересуетесь поддиапазоном интервала Найквиста. Идея позади БПФ изменения масштаба состоит в том, чтобы сохранить то же разрешение, которого вы достигли бы с БПФ в натуральную величину на вашем исходном сигнале путем вычисления маленького БПФ на более коротком сигнале. Более короткий сигнал прибывает из десятикратного уменьшения исходного сигнала. Сбережения прибывают из способности вычислить намного более короткий БПФ при достижении того же разрешения. Это интуитивно: для фактора десятикратного уменьшения D новым уровнем выборки является Fsd = Fs/D, и новым форматом кадра (и длина БПФ) является Ld = L/D, таким образом, разрешение подкошенного сигнала является Fsd/Ld = Fs/L.
DSP System Toolbox предлагает функциональность БПФ изменения масштаба для MATLAB и Simulink через dsp. Системный объект ZoomFFT и блок библиотеки FFT изменения масштаба, соответственно. Следующие разделы обсудят Алгоритм бпф изменения масштаба более подробно.
Прежде, чем обсудить алгоритм используется в 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 и пропустим его через КИХ decimator:
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] сигнала с выборкой Гц Фс уровня. Если мы будем передавать сигнал через комплексный (односторонний) полосовой фильтр, сосредоточенный в ФК = (F1+F2)/2 и с BW пропускной способности = F2 - F1, и затем субдискретизировать его фактором D = пол (Fs/BW), мы понизим желаемую полосу до основной полосы.
В целом, если ФК нельзя выразить в форме k*Fs/D (где K является целым числом), затем переключенный, подкошенный спектр не будет сосредоточен в DC. На самом деле центральная частота ФК будет переведена в [2]:
Fd = ФК - (Fs/D)*floor ((D*Fc + Фс/2) / Фс)
В этом случае мы можем использовать микшер (запускающийся на уровне низкой частоты дискретизации подкошенного сигнала), чтобы сосредоточить желаемую полосу, чтобы обнулить Герц.
Используя пример от предыдущего раздела, мы получаем коэффициенты комплексного полосового фильтра путем модуляции коэффициентов разработанного фильтра 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.
Давайте сначала разработаем одноступенчатый decimator:
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)
Блок FFT изменения масштаба приносит функциональность dsp. ZoomFFT к Simulink. В модели dspzoomfft, мы используем блок FFT изменения масштаба, чтобы осмотреть диапазон частот [800 Гц, 1 600 Гц] входного сигнала, выбранного на уровне 44 100 Гц.
[1] Многоскоростная обработка сигналов - Харрис (Prentice Hall).
[2] Вычисление Переведенных Частот в оцифровке и Субдискретизации Аналоговой Полосы пропускания - Лион (https://www.dsprelated.com/showarticle/523.php)