Другие предметы \ Цифровая обработка сигналов

Методы определения спектральных характеристик электрических сигналов: Учебно методическое пособие, страница 21

Чтобы перейти от индексов к частотной области, надо знать значение шага по времени Ts, через который измерен входной сигнал x(t), и промежуток времени, на протяжении которого он измерялся T – период сигнала. Тогда шаг по частоте в изображении Фурье (т. е. расстояние между соседними гармониками) определяется как:

Df=1/T

%  А диапазон изменения частоты:

F=1/ Ts

% Для рассматриваемого примера: Ts=0.01 с, T=10 с, следовательно, число элементов векторов будет n= 1/Ts * T +1=1001. Вычисляем шаг по частоте и диапазон изменения частоты:

% задаем параметры

Df=0.1;   

F=100;

%  определяем массив частот

f=0: Df: F;

%перестройка вектора частот

f1=- F/2: Df: F/2;

v=fftshift(x);

a=abs(v);

stem(f1,a); grid

title('fourie');

xlabel('frequ, Hz');

%перестройка вектора

частот

f1=- F/2: Df: F/2;

v=fftshift(x);

a=abs(v);

stem(f1,a); grid

title('fourie');

xlabel('frequ, Hz');

Рис.3.Амплитудный спектр сигнала

  

k=1 соответствует нулевое значение частоты, т.е. y(k) содержит значения Фурье-изображения, начиная с нулевой до максимальной частоты F (соответствующей k = n). Функция ejz является периодической по z с периодом 2π. Поэтому информация о Фурье-изображении при отрицательных частотах расположена во второй половине вектора y(k). 

Процедура z=fftshift(y) предназначена для формирования нового вектора z из заданного вектора y путем перестановки второй половины вектора y в первую половину вектора z.

Как и следовало ожидать, огибающая амплитудного спектра представляет собой функцию вида sin(x)/x. Положения нулей огибающей  спектральной плотности определяются длительностью импульса, а расстояние между гармониками – периодом сигнала.

3.  Получение фазового спектра сигнала

                                   Рис.4. Фазовый спектр

fi=angle(v);

stem(f1,fi);grid

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

ПРИЛОЖЕНИЕ 3

Программа  демонстрации эффекта Гиббса

% listing of MAIN.m     % Управляющая программа

clear all

Nf=25; % Number of harmonics % В этой строке задается количество удерживаемых гармоник разложения в ряд Фурье. Меняйте число Nf и наблюдайте, как меняется форма восстановленного сигнала.

k=1:Nf;

T=1;

dt=T/1000;

t=0:dt:T;

A0=AF(0,1)

for k=1 :Nf

    A(k)=AF(k,T);

    B(k)=BF(k,T);

end;

% Short row

Np=1000;

t=0:T/Np:1;

for i=1:Np+1

    S=A0/2;

    for k=1:Nf

        S=S+A(k)*cos(2*pi*k/T*t(i))+B(k)*sin(2*pi*k/T*t(i));

    end;

    s(i)=S;

end;

figure(1), plot(t,s,t,FF(t,T),'r')

% Listing FF.m % Подпрограмма задания сигнала (меандр)

function z=FF(t,T)

% description of signal

N=length(t);

for i=1:N

    if t(i)<0

        z(i)=0;

    end

    if (t(i)>=0)&(t(i)<=T/2)

        z(i)=1/2;

    end;

    if (t(i)>T/2)&(t(i)<=T)

        z(i)=-1/2;

    end;

    if t(i)>T

        z(i)=0;

    end;

end;

% Listing of AF.m % Подпрограмма вычисления  косинусных коэффициентов ряда Фурье

function a=AF(k,T)

% Returns k - coefficients

dt=T/1000;

t=0:dt:T;

F=FF(t,T).*cos(2*pi*k/T*t);

a=2/T*trapz(t,F); % integration

 % Listing of BF.m% Подпрограмма вычисления синусных коэффициентов

     ряда Фурье

function b=BF(k,T)

% Returns k - sin coefficients

dt=T/1000;

t=0:dt:T;

F=FF(t,T).*sin(2*pi*k/T*t);

b=2/T*trapz(t,F); % integration

     На рисунках приведен пример действия программы при восстановлении исходного сигнала по ограниченному числу гармоник Nf  для случаев Nf = 3, 11, 33 и 111.

Скачать файл