PhViewer 2.0

цифровая обработка сигналов

Автор: Степанов А. В., Матвеев С. А.
"Методы компьютерной обработки сигналов систем радиосвязи" 
Все права в отношении данного документа принадлежат автору.

В любой системе электросвязи физическим носителем информации является сигнал, представляющий собой электромагнитное колебание. При этом, управляя при передаче изменением модуляционных параметров сигнала во времени, можно формировать необходимое сообщение. На приемной стороне принятый сигнал может рассматривать как совокупность его мгновенных значений на временном интервале приема сообщения. Известно, что практически любой сигнал s(t) может быть представлен как сумма элементарных колебаний ri(t), умноженных на соответствующим образом подобранные коэффициенты сi 
рис.1
(3.1.1)
Система функций ri(t) носит название базисной системы, а представление сигнала в виде (3.1.1) называют разложением сигнала по сумме базисных функций. Совокупность значений сi называется спектром сигнала в выбранной системе базисных функций. Выбор системы базисных функций для спектрального представления сигналов определяется как объемом априорных знаний об анализируемом сигнале (наличии в нем периодичности, гармоничности и т. д.), так и удобствами практической реализации спектрального анализатора.

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

Математически преобразование Фурье для непрерывной функции х(t) определяется следующим образом 
рис.2
где ω = 2лf, f - частота, 
При этом обратное преобразование Фурье может быть записано в форме 
Спектром дискретного сигнала x(nT) в базисе Фурье называют комплексную функцию 
Очевидно, что расчет спектра по этому выражению требует бесконечного во времени сигнала, что недостижимо при решении практических задач. Обычно имеется только ограниченное число выборок исходного сигнала, по которым и производится оценка его спектра, или, другими словами, осуществляется спектральный анализ. В этом случае под оценкой спектра понимают оценку спектральной плотности сигнала, вкладывая в понятие плотности распределение энергетических и фазовых составляющих ограниченной временной последовательности сигнала по оси частот.

Рассмотрим пример нахождения оценки комплексной спектральной плотности X(jω) дискретного сигнала в произвольной точке f оси частот. Пусть 
где
 - модуль спектральной плотности и φ(ω) =arctg(Ss(ω)/Sc(ω)) - фаза исходного сигнала х(пТ) на угловой частоте ω = 2πf.
Допустим, что сигнал x(nT) ограничен по времени, т. е. он является последовательностью из N отсчетов, взятых с частотой дискретизации fd = 1/Т на интервале от 0 до (N-1)T.

Тогда коэффициенты Ss(ω) и Sc(ω) могут быть найдены из следующих выражений 

Отметим, что существует ряд рекурсивных методов вычисления коэффициентов Ss(ω) и Sc(ω) [18]. Однако наиболее важным при цифровой обработке сигналов является случай, когда необходимо проводить спектральный анализ сигнала не на одной отдельно взятой частоте, а во всей полосе занимаемых им частот. В этом случае результат дискретного преобразования Фурье (ДПФ) сигнала может быть представлен в виде совокупности значений X(fk) в определенных точках fk = k/NT, k = 0,…, N-1 
(3.1.2)

Обратное преобразование для (3.1.2) задается в виде
(3.1.3)

В виду того, что выражение (3.1.2) является линейным преобразованием совокупности значений x(nT), ДПФ имеет вид
(3.1.4) 
где
, WN - матрица размера NхN с элементами WNnk =e-jn2πk/N, n, k = 0,., N-l.
Вычисление ДПФ непосредственно по формуле (3.1.2) требует проведения N2 операций умножения и N(N - 1) операций сложения комплексных чисел. С целью сокращения вычислительных за трат при проведении спектрального анализа разработан ряд алгоритмов, получивших общее название быстрых преобразование Фурье (БПФ). Рассмотрим два базовых алгоритма БПФ для случая, когда N = 2V, V > 0, V — целое.

Алгоритм с прореживанием по времени. Сигнал х(пТ), n= 0,..., N-1 разделяется на два подсигнала с четными и нечетными номерами соответственно, т. е. 

Выражение для ДПФ такого сигнала можно представить в виде
(3.1.5)
где k = 0,..., N-1.
Каждая из сумм в (3.1.5) является результатом ДПФ для N/2 четных и N/2 нечетных отсчетов сигнала х(пТ). Путем следующего аналогичного деления каждого подсигнала на две новых последовательности выражение (3.1.5) может быть представлено виде совокупности преобразований для N/A отсчетов. При последующем делении - для N/8 отсчетов и т. д. пока не останутся только 2-х точечные преобразования. Всего будет v = log2N шагов деления. На каждом m-шаге, m = 0,..., v -1, производится преобразование N комплексных чисел Хm(n) во множество комплексных чисел Хm+1(n) в соответствии со следующим правилом
(3.1.6)
где р, q и r зависят от номера шага m. Эту операцию часто называют «бабочка». 
Графическое представление операции (3.1.6) приведено на рис. 3.1.1. 
Операция «бабочка»
Рис. 3.1.1

Входная последовательность нулевого шага X0n получается перестановкой последовательности х(nT) в соответствии с двоичной инверсией номеров, т. е. х(nT) с номером (nv-1,..., n0) в двоичном представлении запоминается на месте Х0(n0 ,..., nv-1).

Графическая иллюстрация алгоритма БПФ с прореживанием по времени для N = 8 приведена на рис. 3.1.2. 
Алгоритм БПФ с прореживанием по времени для N = 8
Рис. 3.1.2

Пример программы, реализующий алгоритм БПФ с прореживанием по времени, показан ниже.

{число отсчетов N=2v} 
{порядок преобразования PFT=v} 
jj:=1; 
{перестановка} 
for i:=1 to N-1 do 
begin 
if i < jj then 
begin 
x:= fl0][jj]; 
f[0][jj]:= f[0][i]; 
f[0]:= x; 
x:= f[1][jj]; 
f[1][jj]:= f[1][i]; 
f[1][i]:= x; 
end; 
kk:= N div 2; 
while kk < jj do 
begin 
jj:= jj - kk; 
kk:= kk shr 1; 
end; 
jj:= jj + kk; 
end; 
{вычисление коэффициентов} 
for l:= 1 to PFT do 
begin 
s1:= 1 shl l; 
s2:= s1 shr 1; 
x1:= 1; 
yl:= 0; 
z:= cos(Pi/s2); 
t:= -sin(Pi/s2); 
for jj:= l to s2 do 
begin 
i:= jj; 
while i <= N do 
begin 
s0:= I + s2; 
x:= f[0][s0] * xl – f[l][s0] * yl; 
y:= f[l][s0] * x1 + f[0][s0] * yl; 
f[0][s0]:= f[0][i] - x; 
f[1][s0]:= f[1][i] - y; 
f[0][s0]:= f[1][i] + x; 
f[1][i]:= f[1][i] + y; 
i:= i + s1; 
end; 
x:= xl; 
xl:= xl * z – yl * t; 
yl:= yl * z + x * t; 
42 
end; 
end;

В приведенном примере двумерный массив f[j][i] на входе в процесс обработки содержит отсчеты входного сигнала. Он может объявляться и заполняться следующим образом.

Var f: array[0.l, 0.MaxLFT] of single; 
… 
for i:= 0 to N-l do 
begin 
f[0][i]:= InputMasReal[i]; 
f[1][i]:= 0; {= InputMasImage[i]} 
end;

Здесь MaxLFT - максимально ожидаемая длина преобразования, т. е. N ≤ MaxLFT, InputMasReal[i] - массив вещественных отсчетов анализируемого сигнала. Если этот сигнал является комплексным, то вторая половина массива f[1][i] на входе в процесс обработки должна не обнуляться, а заполняться отсчетами мнимой части сигнала, как это показано в комментарии. В результате выполнения БПФ f[0][i] содержит N вещественных коэффициентов (отсчетов) спектральной функции, f[1][i] - N мнимых.

Алгоритм с прореживанием по частоте. Сигнал х(nТ) разбивается на два подсигнала x1(nT) = х(nТ) и х2(nТ) = х((n + N/2)T) так, что 

Далее раздельно вычисляются четные и нечетные отсчеты ДПФ 

Очевидно, что полученные .N/2-точечные последовательности ДПФ аналогичным образом можно представить через N/4-точечные, далее - через N/8-точечные и т. д., пока не останутся только двухточечные. На m-м шаге (m = 0, ..., v-1) производится преобразование N комплексных чисел Xm(n) во множество комплексных чисел Хm+1(n) в соответствии с другой операцией «бабочка», описываемой выражениями 
Xm+1(p) = Xm(p)+X(q), 
Xm+1(q) = Xm(p)-Xm(q)WrN
В алгоритме с прореживанием по частоте, в отличие от алгоритма с прореживанием по времени, в двоично-инверсионном порядке располагается не входная последовательность временных отсчетов сигнала х(nT), а выходная последовательность значений спектральных составляющих Х(k). Графическая иллюстрация алгоритма БПФ с прореживанием по времени для N= 8 приведена на рис. 3.1.3. 
Алгоритм БПФ с прореживанием по времени для N= 8
Рис 3.1.3

Приведенные алгоритмы БПФ предпочтительнее по числу выполняемых операций примерно в N/v раз по сравнению с прямым методом вычисления ДПФ по формуле (3.1.2). Существует целый ряд других алгоритмов БПФ, оптимизированных по объему вычислений и требуемой памяти для различных значений N. Их подробное описание с анализом вычислительной сложности можно найти в [19, 20].

При практической реализации различных алгоритмов БПФ необходимо учитывать следующие существенные факторы. Если объем доступной для анализа последовательности отсчетов сигнала меньше требуемого по условиям выполнения БПФ значения N, то ее можно дополнить в конце нулевыми отсчетами. Частотный диапазон анализа в общем случае определяется частотой Найквиста (см. раздел 2.2) и составляет 0.1/2Т Гц. В этом случае предельная разрешающая способность анализа по частоте Δf задается следующим выражением 
Δf = m/NT, (3.1.7)

где величина m ∈ [0.1] определяется параметрами дополнительной обработки исходной временной и (или) полученной частотной последовательности отсчетов сигнала.
Точное равенство m = 1 в (3.1.7) достижимо только при отсутствии такой обработки, но в этом случае значения реально вычисленной спектральной функции, как правило, будут отличаться от теоретически ожидаемых значений. Рассмотрим причины такого несовпадения.

Конечная длина выборки входного сигнала в сочетании с его аппроксимацией при ДПФ линейной комбинацией функций синусов и косинусов с кратными частотами вызывает эффект размывания спектральной функции. Иначе говоря, этот эффект обусловливает концентрацию спектральной функции не только в значимых точках частотной оси, но и распределение функции в смежных с ними областях. Графически этот эффект для случаев «идеально» и реально вычисленных спектров синусоидального сигнала показан на рис. 3.1.4. 
Эффект размывания спектральной функции
Рис. 3.1.4

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

Сглаживание временной последовательности заключается в умножении отсчетов сигнала на весовые коэффициенты специальной функции, называемой «окном». Коэффициенты оконной функции подбираются с целью уменьшения влияния отсчетов на концах последовательности на результат вычисления ДПФ. Фактически это является попыткой преодоления ограниченности объема обрабатываемой выборки. Разработано большое количество оконных функций, отличающихся по совокупности параметров разрешающей способности, степени сглаживания, ухудшения отношения сигнал/шум после применения оконной функции и т. д. Не вдаваясь в теоретические характеристики и особенности применения отдельных окон, описание которых можно найти, например, в [21], приведем только наиболее часто используемые при решении практических задач оконные функции pwn(t), упорядоченные по степени сглаживания. 
1. Прямоугольное окно. Для этого окна 
pwN[nT] = 1, n = 0…N-1.
Практически такое окно эквивалентно отсутствию дополнительной обработки временной последовательности, что приводит к m = 1 в (3.1.7). Окно целесообразно использовать в случае, когда гармонические составляющие анализируемой последовательности совпадают с сеткой частот ДПФ (см. (3.1.2)).

2. Косинусное окно 
Косинусное окно

где [.] — операция взятия целой части числа и N > 10.
Это окно сглаживает крайние отрезки временного ряда, составляющие примерно по 10% от его длины. В (3.1.7) для такого окна m ≈ 0,9.

3. Обобщенное окно Хэмминга 
pwN[nT] = α - (1-α)cos(2πnT),   α∈]0...1[, n = 0...N-1
Обычно α = 0,54. Для такого случая в (3.1.7) m ≈ 0,55.

Кроме оконных функций, применяемых во временной области, существуют различные окна, сглаживающие результаты расчета спектральной функции. Простейшим таким окном является сглаженная оценка спектра, полученная усреднением L соседних отсчетов. Тогда один отсчет энергетического спектра сигнала , где f = 1/NT, k∈[0…N/2], можно оценить как 
Полученное выражение можно заменить сверткой более общего вида 

В этом случае часто используют биноминальные сглаживающие веса а(l), которые вычисляются по формуле а(l) = ã(l)/2L-1 где ã(l) принимает значения (1,2,1) при L = 3, (1,3,3,1) при L = 4, (1,4,6,4,1) при L = 5 и т. д.

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

{предварительный расчет коэффициентов оконной функции} 
for i:= 0 to N-l do 
pw[i]:= 0.54 - 0.46 * cos(i*2*Pi/N); 
… 
{сглаживание входных данных} 
for i:= 0 to N-l do 
begin 
fl0][i]:= InputMasReal[i] * pw[i]; 
f[1][i]:= 0; 
end; 
… 
{расчет энергетического спектра} 
for i:= 0 to N div 2 - 1 do 
aa[i]:= sqrt(f[0][i] * f[0][i] + f[l][i] * f[1][i]); 
{сглаживание энергетического спектра} 
m = L div 2; 
for k:= 0 to m do 
begin 
j:= N div 2 - 1 - k; 
for i:= 0 to j do aa[i]::= aa[i] + aa[i + 1]; 
for i:= 0 to j-1 do aa[j - i]:= aa[j - i] + aa[j – i - l]; 
aa[0]:= 2*aa[0] * aa[0]; 
end;

В качестве альтернативы проведению спектрального анализа дискретных гармонических сигналов на основе ДПФ можно использовать спектральный анализ на основе дискретного преобразования Хартли [22]. Кроме того, существуют специальные системы базисных функций для других типов дискретных сигналов, например бинарных [23, 24]. Также возможно проведение спектрального анализа другими методами: фильтрационным, корреляционным и т. д. Не вдаваясь в подробности их описания, еще раз отметим, что основной целью спектрального анализа, независимо от способа его реализации, является получение спектральной плотности сигнала X(jω) = S(ω)еjφ(ω). Наличие оценки спектральной плотности позволяет, в свою очередь, реализовывать различные процедуры обнаружения сигнала, анализа его технических параметров, помехоустойчивой обработки и т. д.

В качестве примера одной из задач, возникающих в процессе спектрального анализа, рассмотрим задачу определения вида модуляции сигнала.

Пусть на входе спектрального анализатора может присутствовать один из сигналов следующих типов: двухпозиционный ЧM-сигнал с индексом манипуляции больше 1, двухпозиционный ЧМ-сигнал с индексом манипуляции равным 0,5, т. е. ММС-сигнал, а также двух- и четырехпозиционный ФМ-сигналы. Усредненные энергетические спектры этих сигналов показаны соответственно на рис. 3.1.5, а, б, в, г.
Усредненные энергетические спектры
Рис. 3.1.5

Их визуальный анализ показывает, что только сигнал ЧМ2 может однозначно определяться по своему спектру. Идентификацию остальных сигналов можно провести, используя операцию умножения фазы, часто применяемую при демодуляции фазовых сигналов.

Известно, что умножение текущей фазы ФМ-сигнала на 2, 4 и т. д. позволяет последовательно снимать по одному уровню манипуляции. Простым способом такого умножения является возведение сигнала в соответствующую степень. Пусть, например, анализируемый сигнал описывается следующим выражением s(t) =Acos(ωt + θ), где в случае обработки ФМ2-сигнала θ = {0, π}. Возведем этот сигнал в квадрат 
s2(t) = A2cos2 (ωω+ θ) = A2(1/ 2 +1/ 2cos( 2ωt + 2θ)) = A2(1/ 2 +1/ 2cos( 2ωt)) .
Таким образом, мы получили выражение, описывающее сигнал с удвоенной частотой, и не содержащее информацию об исходной фазе сигнала. Очевидно, подобный подход можно использовать для снятия манипуляции у ФМ4, ФМ8 и т. п. сигналов.

Каждая операция возведения в степень приводит к соответствующей трансформации спектра анализируемого сигнала. В нем возникают гармоники удвоенной, учетверенной и т. д. центральной частоты, которые и являются признаками, характеризующими конкретный тип сигнала. На рис. 3.1.6 и 3.1.7 показаны энергетические спектры рассматриваемых сигналов, полученные после возведения этих сигналов во вторую и соответственно четвертую степени.
Рис. 3.1.6 
Рис. 3.1.7

Визуальный анализ приведенных спектров показывает, что, как и ожидалось, для сигнала ФМ2 гармоника удвоенной частоты появляется после возведения сигнала в квадрат, а для ФМ4-сигнала гармоника учетверенной центральной частоты сигнала возникает после возведения сигнала в четвертую степень. Возведение в квадрат сигнала ММС приводит к появлению двух гармоник, расположенных симметрично относительно удвоенной центральной частоты сигнала, с разносом между ними, равным удвоенному исходному разносу.

Программирование операции возведения сигнала в соответствующую степень является очевидным. Такую процедуру можно реализовать, например, следующим образом

{возведение в степень и сглаживание входных данных 
for i:= 0 to N-l do 
begin 
f[0][i]:= lntPower(InputMasReal[i], iSpectraPower) * pw[i]; 
f[1][i]:= 0; 
end;

Здесь переменная iSpectraPower является требуемой степенью.

Завершая рассмотрение методов спектрального анализа, следует сказать следующее. Спектральный анализ благодаря относительно просто реализуемой возможности комплексного представления сигналов во временной и частотных областях является, по-видимому, одним из наиболее часто используемых инструментов в теории и практике цифровой обработки сигналов. В следующих разделах будет показаны возможности его применения не только для непосредственного анализа параметров сигнала, но и для реализации квадратурной формы его представления, при разработке цифровых фильтров с постоянными и адаптивно подстраиваемыми параметрами, при помехоустойчивой обработке речевых сигналов и т. д.