Регулярный БПФ

Разрешение сигнала ограничено его длиной. Мы проиллюстрируем этот факт простым примером. Считайте сигнал сформированным суммой двух синусоид:

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 Системный объект и библиотечный блок БПФ изменения масштаба, соответственно. Следующие разделы обсудят Алгоритм бпф изменения масштаба более подробно.

Подход микшера

Прежде, чем обсудить алгоритм использовало в dsp.FFT, мы представляем подход микшера, который является популярным методом БПФ изменения масштаба.

Для примера здесь, примите, что вы только интересуетесь интервалом [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 спроектирует фильтр и применит его к входному сигналу.

Давайте использовать 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)

Масштабируйте блок Simulink БПФ

Блок FFT изменения масштаба приносит функциональность dsp.ZoomFFT к Simulink. В модели dspzoomfft, мы используем блок FFT изменения масштаба, чтобы смотреть диапазон частот [800 Гц, 1 600 Гц] входного сигнала, произведенного на уровне 44 100 Гц.

Ссылки

[1] Многоскоростная обработка сигналов - Харрис (Prentice Hall).

[2] Вычисление Переведенных Частот в оцифровке и Субдискретизации Аналоговой Полосы пропускания - Лион (https://www.dsprelated.com/showarticle/523.php)