Алгоритм БПФ предназначен для одновременного расчёта всех спектральных отсчётов X(n). Если же необходимо получить эти отсчёты лишь для некоторых n, может оказаться предпочтительнее применение алгоритма Герцеля.
Функция передачи отдельного канала ДПФ [3, c.262]:
,
(2)
где n – номер текущего (анализируемого) отсчёта, N – количество отсчётов в выборке.
Эта форма соответствует рекурсивному фильтру N-го порядка. Поэтому вычисление отдельного спектрального отсчёта можно реализовать с помощью такого фильтра.
Однако, для удобства практического использования формулу (2) можно модифицировать. Прежде всего, заметим, что слагаемое –z-N в числителе выполняет лишь роль «ограничителя», заставляя импульсную характеристику закончиться, достигнув длины N отсчётов. Такой фильтр вычисляет отсчёт мгновенного спектра сигнала, всё время, используя последние N отсчётов входного сигнала. Если же мы будем обрабатывать сигнал порциями по N отсчётов, перед каждой порцией обнуляя состояние фильтра, то конечность длины импульсной характеристики не будет иметь значения. В этом случае можно удалить слагаемое –z-N из числителя функции передачи, получив, таким образом, рекурсивный фильтр первого порядка:
(3)
Следующее изменение связано с тем, что коэффициент обратной связи в единственной рекурсивной ветви данного фильтра является комплексным. Поэтому при обработке каждого входного отсчёта необходимо выполнить четыре вещественных умножения и столько же вещественных сложений. Число операций можно сократить, умножив числитель и знаменатель функции передачи на (1-exp(-j·2π·n/N)· z-1), перейдя за счёт этого к фильтру второго порядка (4).
.(4)
Рекурсивная часть получившегося фильтра имеет вещественные коэффициенты, к тому же коэффициент при z-2 равен единице. Поэтому при обработке каждого входного отсчёта в рекурсивной ветви фильтра необходимо выполнить всего одно вещественное умножение и два вещественных сложения.
В рекурсивной ветви фильтра комплексные коэффициенты исчезли, но они, естественно, появились в числителе функции передачи, то есть в нерекурсивной части фильтра. Однако, если реализовать фильтр в канонической форме, то расчёты будут выполняться сначала в рекурсивной ветви, а потом в нерекурсивной. Полученная структура фильтра показа на рисунке 6. Данный способ вычисления спектрального отсчёта называется алгоритмом Герцеля.
При разложении функции в ряд Фурье полагаем, что функция периодическая, с периодом, равным N·Δt, где Δt – период выборок. При этом на границах периодов такая функция будет иметь разрывы, так как исходная функция не была периодической с периодом длительности последовательности выборок. Разрывы в функции сильно отражаются на ее спектре, искажая его. Для устранения этого эффекта применяются так называемые взвешивающие окна. Они плавно снижают амплитуду сигнала вблизи краев анализируемого участка. Весовые окна имеют форму, похожую на гауссиан. Выбранный для анализа участок сигнала домножается на весовое окно, которое устраняет разрывы функции при «зацикливании» данного участка сигнала. «Зацикливание» происходит при ДПФ, так как алгоритм ДПФ полагает, что функция периодическая.
Важное свойство спектрального анализа заключается в том, что не существует одного, единственно правильного спектра какого-либо сигнала. Спектр можно вычислять с применением различных размеров ДПФ и различных весовых окон. Для каждого конкретного приложения предпочтительно использовать свои способы.
Еще одно важное свойство заключается в том, что при разложении в спектр мы находим не те синусоидальные составляющие, из которых состоял исходный сигнал, а лишь находим, с какими амплитудами нужно взять определенные кратные частоты, чтобы получить исходный сигнал.
Рисунок 6 – Алгоритм Герцеля
Другими словами, разложение проводится не по «частотам источника», а по «частотам алгоритма БПФ». Однако обычно (особенно при использовании весовых окон) этого почти не заметно по графику спектра, то есть график спектра достаточно адекватно отображает именно частоты исходного сигнала.
2 Реализация алгоритма Герцеля
Для реализации алгоритма Герцеля разобьём фильтр (рисунок 6) на рекурсивную и трансверсальную части (рисунок 7).
Рисунок 7 – Алгоритм Герцеля (каскадная реализация): рекурсивная и трансверсальная части
Рекурсивная часть должна выполняться N раз.
Рисунок 8 – Алгоритм Герцеля: получение одного спектрального отсчёта
На рисунке 8 Ftone – искомый спектральный отсчёт (частота, Гц); FS – частота получения выборок (сэмплирования), Гц; N – количество элементов в выборке; n – номер текущего элемента; input – очередной входной элемент выборки (с номером n); q1 – первый элемент задержки; q2 – второй элемент задержки; q0 – текущее значение на выходе рекурсивной ветви.
Поскольку нас интересует только конечный результат вычисления спектрального отсчёта, комплексное умножение и сложение в нерекурсивной ветви должны быть выполнены только один раз – для получения окончательного результата.
Так как нас интересует только амплитудный спектр, и не интересует фазовый, то мы после нахождения комплексного выходного отсчёта в трансверсальной части находим корень из квадратов мнимой и реальной частей нужного спектрального отсчёта (амплитуду). Алгоритм нахождения одного спектрального отсчёта изображён на рисунке 8 в виде блок-схемы.
В случае необходимости получения нескольких спектральных отсчётов с определённым шагом потребуются элементы задержки в виде массивов, размерность которых будет соответствовать количеству получаемых спектральных отсчётов. Массивы необходимо обнулять до получения первого отсчёта.
Например, функция, реализующая вычисление двухсот отсчётов амплитудно-частотного спектра (АЧС) от 0 до 398 Гц с шагом 2 Гц, записанная на языке Си, может выглядеть следующим образом:
#include <math.h>
#define SAMPLING_RATE 100 // частота получения выборок
#define STEP 2 // шаг (расстояние между спектральными отсчётами)
#define MAX_BINS 200 // количество вычисляемых спектральных отсчётов
#define GOERTZEL_N 4096 // количество точек в выборке, на базе которой вычисляется АЧС
#define pi 3.141592654
double retZero = -1.0;
unsigned int sample_count=0;
double* goertzel( unsigned int sample )
{
double q1[ MAX_BINS ]={0.0};
double q2[ MAX_BINS ] ={0.0};
double r[ MAX_BINS ] ={0.0};
double q0;
unsigned int i;
sample_count++;
for ( i=0; i<MAX_BINS; i++ )
{
q0 = 2.0 * cos(2.0 * pi * ((double)(i* STEP))/SAMPLING_RATE) * q1[i] - q2[i] + sample;
q2[i] = q1[i];
q1[i] = q0;
}
if (sample_count == GOERTZEL_N)
{
for ( i=0; i<MAX_BINS; i++ )
{
r[i] = (q1[i] * q1[i]) + (q2[i] * q2[i])
- (2.0 * cos(2.0 * pi * ((double)(i* STEP))/SAMPLING_RATE) * q1[i] * q2[i]);
r[i] = sqrt(r[i]);
q1[i] = 0.0;
q2[i] = 0.0;
}
sample_count = 0;
return r;
}
return &retZero;
}
Как только функция вернёт указатель на положительное число – то это будет являться знаком, что получена очередная порция спектральных отсчётов АЧС. Изменяя константу MAX_BINS можно получить требуемое число точек, например, для отображения АЧС на жидкокристаллическом индикаторе (ЖКИ).
Дальнейшей разработкой является универсализация функции для получения спектральных отсчётов АЧС, рассчитанных для произвольных начальной частоты, шага, частоты дискретизации, и количества точек в выборке, которое, в данном случае, не обязательно может представляться как 2m, а быть произвольным целым числом.
Ниже приведён код усовершенствованной реализации алгоритма Герцеля с упомянутыми выше свойствами.
#include <math.h>
#define PI 3.141592654
#define MAX_BINS 200
double* goertzel(int* sample, float Fstart ,double step, unsigned long Fsampling, unsigned int GoertzN)
{
double q1[ MAX_BINS ] ={0.0};
double q2[ MAX_BINS ] ={0.0};
double r[ MAX_BINS ] ={0.0};
double q0;
unsigned int i;
for(unsigned int k=0; k<GoertzN; k++){
for ( i=0; i<MAX_BINS; i++ )
{
q0 = (2.0 * cos(2.0 * PI * ((double)Fstart+i*((double)step))/(double)Fsampling) * q1[i]) - q2[i] +* sample;
q2[i] = q1[i];
q1[i] = q0;
}
sample++;
if (k == GoertzN-1)
{
for ( i=0; i<MAX_BINS; i++ )
{
r[i] = sqrt((q1[i] * q1[i]) + (q2[i] * q2[i]) - (2.0 * cos(2.0 * PI * ((double)Fstart+i*((double)step))/(double)Fsampling) * q1[i] * q2[i]));
q1[i] = 0.0;
q2[i] = 0.0;
}
sample_count = 0;
}
}
return r;
}
На этом этапе дальнейшая разработка сопряжена с тестированием алгоритма в работе и дальнейшим его усовершенствованием путём ускорения вычислений, устранения нежелательных эффектов (например, боковых лепестков из-за конечности выборки и непериодичности сигнала). Этому посвящены следующие пункты, относящиеся к разработке методики анализа спектра, основываясь на алгоритме Герцеля.
3 Тестирование алгоритма Герцеля
Для удобного тестирования полученного алгоритма автором была написана специальная программная оболочка на WinAPI, пользовательский графический интерфейс которой представлен на рисунке 9.
К программе на этапе препроцессорной обработки подсоединяется модуль с функцией, реализующей алгоритм получения спектральных отсчётов. В качестве входных параметров берутся указанные пользователем данные в текстовых полях: диапазон просматриваемых частот, количество отсчётов, считываемых из файла сигнала для анализа и частота семплирования. На выходе программа генерирует m-файл, который содержит вектор со спектральными отсчётами, полученными в результате исполнения алгоритма. Полученный m-файл легко загружается в MATLAB для анализа, например для построения на основе его графиков.
Рисунок 9 –
Интерфейс программы тестирования алгоритма Герцеля
В дальнейшем программа была дополнена процедурой измерения времени исполнения фрагмента кода. Для более точного измерения функция вычисления спектра запускается 100 раз, затем проводится усреднение времени однократного выполнения:
clock_t time1,time2;
time1 = clock();
for(int zzz=0; zzz<100; zzz++)
OutDouble = goertzel(SignalInput, Fstart, step, Fsampling, N_numeric);
time2 = clock();
char strValueMilliseconds[18];
sprintf(strValueMilliseconds,"%ld.%d мс",((time2-time1)*10.0/CLOCKS_PER_SEC),
((time2-time1)*10.0/CLOCKS_PER_SEC)%100);
MessageBox(hwndMain, strValueMilliseconds,"Working time is",
MB_OK|MB_ICONINFORMATION|MB_APPLMODAL);
Данное дополнение потребуется для относительной оценки изменений во времени исполнения алгоритма в зависимости от производимых изменений в коде с целью получить лучший по качеству и быстродействию код.
Рисунок 10 – Диалоговое окно оценки времени выполнения алгоритма
При этом тестирование на относительное быстродействие проводилось на одной и той же конфигурации – использовался специальный отдельный персональный компьютер на базе процессора Intel Atom 230 и чипсета Intel 945GC.
Для тестирования алгоритма этой программой с помощью MATLAB было сгенерировано 60 специальных файлов. Каждый из этих файлов представляет эмулированные отсчёты 12-и разрядного АЦП для синусоидального сигнала различных частот:
script
fd=100; % условная частота дискретизации
for j= 1:1:60
for i= 1:1:4096
X(i,j)=2048+fix(2047*sin((2*pi/100)*i*j/2)); % эмулируем отсчеты АЦП
end
end
dlmwrite('file001.sig', X(:,1)); % Пишем файлы
dlmwrite('file002.sig', X(:,2));
dlmwrite('file003.sig', X(:,3));
dlmwrite('file004.sig', X(:,4));
dlmwrite('file005.sig', X(:,5));
dlmwrite('file006.sig', X(:,6));
dlmwrite('file007.sig', X(:,7));
dlmwrite('file008.sig', X(:,8));
…
dlmwrite('file060.sig', X(:,60));
Зададимся входными параметрами (рисунок 9):
- диапазон просматриваемых частот: от 0 до 30 Гц;
- длина выборки сигнала 4096 точек;
- частота семплирования 100 Гц;
- входной файл сигнала file020.sig (частота синусоидального сигнала 10 Гц (10 точек
на период)).
В результате исполнения программы был получен файл out1.m:
script
Yout=[
5095.000000;
4410.679465;
2614.523392;
…
1662.687095];
На основании данных из полученного файла строим нужный нам график АЧС:
>> out1
>> Z=20*log10(Yout);
>> x=0:30/199:30;
>> plot(x,Z-max(Z));
>> grid on;
График представлен на рисунке 11.
При этом наблюдаются значительные боковые лепестки на уровне до -13 дБ.
Выполнив аналогичный расчёт для другого сигнала: file022.sig, частота 11 Гц, получим спектральную картину, изображённую на рисунке 12.
Как видим, уровень боковых лепестков полученного АЧС для этого сигнала значительно ниже, и составляет порядка -36 дБ, т.е. спектр более «чистый».
Причина появления больших боковых лепестков кроится не в несоответствии периодичности сигнала длительности выборки. Если посмотреть на начальные и конечные участки первого и второго сигналов, можно заметить, что синусоиды усечены на разных уровнях (смотрите рисунок 13 и рисунок 14).
Рисунок 11 – График АЧС 200-точечного ДПФ методом Герцеля для файла file020.sig
Однако при этом разница начальной и конечной фаз в обоих случаях отличаются не сильно.
На самом деле, причина образования боковых лепестков здесь в другом. При разложении в спектр мы находим не те синусоидальные составляющие, из которых состоял исходный сигнал, а лишь находим, с какими амплитудами нужно взять определенные кратные частоты, чтобы получить исходный сигнал, раскладывая его, по сути, по частотам БПФ.
Применяя этот же алгоритм, но с шагом по частоте, соответствующим шагу для БПФ, мы получим участок спектральной картинки, практически соответствующий той, что была получена при применении классического БПФ.
Рисунок 12 – График АЧС 200-точечного ДПФ методом Герцеля для файла file022.sig
Рисунок 13 – Начальный и конечный участки выборки для сигнала из file020.sig
Рисунок 14 – Начальный и конечный участки выборки для сигнала из file022.sig
Шаг по частоте в спектральном анализе БПФ рассчитывается по формуле:
(5)
Полученная с помощью алгоритма Герцеля спектральная картинка, при применении начальной частоты вычисления кратной ΔF и шага по частоте равному ΔF для сигнала из файла file020.sig представлена на рисунке 15.
Как видим, никаких побочных лепестков в АЧС не наблюдается. Вычисления с помощью этого алгоритма и построения графиков были проведены автором для всех шестидесяти сигналов, значительных боковых лепестков при этом обнаружено не было.
Для уменьшения боковых лепестков в случае вычислений АЧС с произвольным начальным значением и шагом по частоте и уменьшения эффекта растекания спектра в случае использования по частоте равному или кратному FS/N можно применить специальное временное окно.
Ограничивая выборку определённым количеством точек, мы вызываем растекание спектра (spectral leakage), этим вызвано то, что на АЧС мы видим не спектральную линию, соответствующую частоте синусоидального сигнала, а нарастающую до максимума и спадающую линию как на рисунке 15. То есть мы непроизвольно используем прямоугольное окно (называемое в научной литературе также окном Дирихле). И чем короче наша выборка (меньшее количество точек), тем сильнее эффект растекания спектра.
Рисунок 15 – График АЧС ДПФ методом Герцеля для файла file020.sig при шаге по частоте, равном FS/N
Растекание спектра от сильных сигналов может скрыть под собой слабые сигналы, находящиеся близко по частоте. Прямоугольное окно даёт хорошую разрешающую способность по частоте для сигналов, не сильно различающихся по мощности, но не для сигналов с сильно различающимися амплитудами. Т.о. для сигналов, близких по частоте будет иметь место малый динамический диапазон (ДД) по амплитуде.
Применение оконных временных функций даёт расширение ДД за счёт уменьшения эффекта растекания спектра, однако, также даёт и некоторое снижение чувствительности и расширение главного лепестка АЧС.
Существует множество весовых окон, отличных от прямоугольного, названных в честь их создателей.
Как выяснилось в результате тестирования, окно Блэкмана наиболее подходит для снижения уровня боковых лепестков в данном случае.
Окно Блэкмана по форме напоминает гауссову функцию (рисунок 16):
>> x=blackman(1000000);
>> plot(x);
>> grid on;
Рисунок 16 – Окно Блэкмана
Классическое выражение для вычисления имеет форму (6), эта же форма использована и в среде MATLAB:
(6)
Свертывая временные отсчёты сигнала с оконной функцией, мы теряем в среднем его энергетику. Выясним, на сколько требуется предварительно умножить все временные отсчёты перед применением временного окна Блэкмана:
>> y=0;
>> for i=1:30000000
y=y+x(i);
end
>> 30000000/y
ans =
2.38095238095238
Таким образом, итоговое выражение для функции окна будет иметь вид:
(7)
Функция для вычисления отсчётов для временного окна может быть представлена как:
double blackman(int * sample, const unsigned int GoertzN, const unsigned int i ){
double z=(float)(*sample)-2048.;
z *= 1.-1.19047619047619 *cos(2*PI*i/(GoertzN-1))+ 0.19047619047619*cos(4*PI*i/(GoertzN-1));
return z;
}
Здесь аргументы функции blackman(): указатель на значение очередного входного отсчёта сигнала sample, количество отсчётов выборки GoertzN, номер текущего отсчёта i.
Соответствующая часть алгоритма Герцеля соответственно изменится:
for(unsigned int k=0; k<GoertzN; k++){
for ( i=0; i<MAX_BINS; i++ )
{
q0 = (2.0 * cos(2.0 * PI * ((double)Fstart+i*((double)step))/(double)Fsampling) *
q1[i]) - q2[i]
+ blackman(sample, GoertzN, k);
q2[i] = q1[i];
q1[i] = q0;
}
…
Полученная с помощью алгоритма Герцеля спектральная картинка с применением окна Блэкмана при применении начальной частоты вычисления кратной ΔF и шаге по частоте, равном FS/N для сигнала из файла file020.sig представлена на рисунке 17.
Рисунок 17 – График АЧС ДПФ методом Герцеля
для файла file020.sig
при шаге по частоте, равном FS/N
и применении к временным отсчётом окна Блэкмана
Быстродействие алгоритма также снизилось и теперь составляет в среднем 587,98 мс. против 174,37 мс. без окна.
Применение окна для случая с произвольным шагом по частоте также показывает хорошие результаты. Например, при тех же входных данных, для которых был получен график АЧС, изображенный на рисунке 11, с применением окна Блэкмана получим спектральную картинку, изображённую на рисунке 18.
Рисунок 18 – График АЧС при ДПФ методом Герцеля для файла file020.sig при шаге по частоте, не кратном FS/N и применении к временным отсчётом окна Блэкмана
Как видим, уровень боковых лепестков резко снизился.
4 Повышение быстродействия алгоритма
Алгоритм в разработанном на данный момент виде позволяет вычислять АЧС для нужного количества точек с произвольным шагом. При этом сохраняется исходный массив выборок и расход дополнительной памяти весьма мал. Однако, введение временного окна значительно увеличило время вычислений, поскольку применение его не оптимизировано и его функция вычисляется каждый раз, когда для вычислений требуется считать значение временного отсчёта.
Введя дополнительную переменную, можно снизить трудоёмкость алгоритма. Введем, например переменную windowedSample:
…
double windowedSample;
for(unsigned int k=0; k<GoertzN; k++){
windowedSample=blackman(sample,GoertzN,k);
for ( i=0; i<MAX_BINS; i++ )
{
q0 = (2.0 * cos(2.0 * PI * ((double)Fstart+i*((double)step))/(double)Fsampling) * q1[i]) - q2[i]
+windowedSample;
q2[i] = q1[i];
q1[i] = q0;
}
sample++;
…
В результате быстродействие значительно возросло: среднее время выполнения стало 236,54 мс.
Очевидно, что для ещё большего увеличения быстродействия можно, если позволяет запас по ПЗУ заранее вычислить все нужные коэффициенты окна. Причём, т.к. окно симметричное, то достаточно хранить половину его коэффициентов. Правда, это накладывает ограничения на количество точек исходной выборки, для снятия которого может потребоваться применение интерполяции. Однако, в большинстве случаев представляется возможность выбора количества выборок сигнала как 2n (например, 1024, 2048, 4096…)
Вычислим коэффициенты окна Блэкмана для числа точек N = 4096:
>> X=blackman(4096);
Z=4096/sum(X);
X=X*Z;
dlmwrite('blackman.c', X, 'precision', ‘%e’);
Запишем половину получившегося окна в виде массива на языке Си:
const double BlackmanWindow[2048]={
-3.305042087377649e-017,
5.046050058672026e-007,
2.018425171699819e-006,
4.541475940505306e-006,
8.073783051521783e-006,
1.261538253803671e-005,
1.816632072884110e-005,
2.472665424591600e-005,
3.229645000439942e-005,
………….
2.381404508279874e+000,
2.381436689488475e+000,
2.381464273686991e+000,
2.381487260734456e+000,
2.381505650513398e+000,
2.381519442929837e+000,
2.381528637913288e+000,
2.381533235416761e+000};
В результате функция вычисления значений оконной функции значительно упрощается:
double blackman(int * sample, unsigned int GoertzN, unsigned int i ){
return (i<2048)? *sample* BlackmanWindow[i]: *sample* BlackmanWindow[4095-i];
}
Среднее время вычисления двухсот точек АЧС снизилось до 210,78 мс.