Регулярное БПФ

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

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

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

Блок БПФ приносит функциональность dsp. Масштабирование FFT до Simulink. В модели dspzoomfft мы используем блок БПФ для просмотра полосы частот [800 Гц, 1600 Гц] входного сигнала, дискретизированного с частотой 44100 Гц.

Ссылки

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

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