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

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

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, а новый размер кадра (и длина FFT) Ld = L/D, поэтому разрешение прореженного сигнала Fsd/Ld = Fs/L.

Панель инструментов DSP System Toolbox предлагает функции масштабирования FFT для MATLAB и Simulink через dsp. Объект ZoomFFT System и блок библиотеки масштабирования FFT соответственно. В следующих разделах более подробно рассматривается алгоритм масштабирования БПФ.

Подход смесителя

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

В данном примере предположим, что вас интересует только интервал [16 Гц, 48 Гц].

BWOfInterest = 48 - 16;
Fc           = (16 + 48)/2; % Center frequency of bandwidth of interest

Достижимый коэффициент прореживания:

BWFactor = floor(Fs/BWOfInterest);

Подход микшера состоит в том, чтобы сначала сдвинуть интересующую полосу до постоянного тока с помощью микшера, а затем выполнить фильтрацию нижних частот и прореживание с коэффициентом BWFactor (с использованием эффективной полифазной структуры децимации FIR).

Давайте разработаем коэффициенты прореживающего фильтра с помощью функции designMultirateFIR:

B = designMultirateFIR(1,BWFactor);
D = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',B);

Теперь давайте смешаем сигнал до постоянного тока и фильтруем его через дециматор FIR:

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

Теперь возьмем БПФ отфильтрованного сигнала (заметим, что длина БПФ уменьшается на БПФ, или длина прореживания, по сравнению с обычной БПФ, при сохранении того же разрешения):

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)

В этом случае можно использовать смеситель (работающий при низкой частоте дискретизации прореженного сигнала) для центрирования желаемой полосы до нуля Герца.

На примере из предыдущего раздела получаем коэффициенты комплексного полосового фильтра, модулируя коэффициенты проектируемого фильтра нижних частот:

Bbp     = B .*exp(1j*(Fc / Fs)*2*pi*(0:numel(B)-1));
D = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',Bbp);

Теперь выполним фильтрацию и FFT:

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] Гц. Достижимым коэффициентом прореживания является floor (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 является трехступенчатым многоскоростным фильтром нижних частот. Преобразуем его в полосовой фильтр, выполняя сдвиг частоты коэффициентов каждого этапа, принимая во внимание коэффициент кумулятивной прореживания:

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 - это системный объект, реализующий масштабирование FFT на основе многоступенчатого многоступенчатого полосного фильтра, подсвеченного в предыдущем разделе. Задаются требуемая центральная частота и коэффициент прореживания, а также 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 мы используем блок масштабирования БПФ для проверки полосы частот [800 Гц, 1600 Гц] входного сигнала, дискретизированного на частоте 44100 Гц.

Ссылки

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

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