ЦИФРОВАЯ ОБРАБОТКА СИГНАЛОВ
Работа добавлена: 2016-05-20





Минобрнауки России

Федеральное государственное бюджетное общеобразовательное

учреждение высшего профессионального образования

«Тульский государственный Университет»

Институт высокоточных систем им. В. П. Грязева

КАФЕДРА РАДИОЭЛЕКТРОНИКИ

ПОЯСНИТЕЛЬНАЯ ЗАПИСКА

к контрольно-курсовой работе по дисциплине

«ЦИФРОВАЯ ОБРАБОТКА СИГНАЛОВ».

 Выполнил студент гр. 130611 _________________ Нгуен Суан Чыонг.

 Проверил: Доцент                ___________________ Полынкин А. В.

Тула 201 5

Аннотация

В даннойконтрольно-курсовой работе содержат 6 заданий в процессе выполнении лабораторной работы по предметам ЦОС. Произведены результаты исследования, вычисления заданий на самостоятельную работу. Расчёты проводятся с помощью компьютера с использованием программы “MATLAB 2013a”.

Пояснительная записка содержит 53 листов и 52 рисуноков.

Содержание

  • Аннотация
  • Задание 1: Дискретное преобразование Фурье (часть 1)
  • Задание 2: Дискретное преобразование Фурье (часть 2)
  • Задание 3: Спектральный анализ: непараметрические методы
  • Задание 4: Спектральный анализ: параметрические методы
  • Задание 5: Многоскоростные системы ЦОС
  • Задание 6: Адаптивные фильтры
  • Заключение
  • Список использованной литературы

Задание 1: Дискретное преобразование Фурье (часть 1)

Исходные амплитуды и частоты дискретных гармоник:

A1 = 1.015   A2 = 2.03   f1 = 250   f2 = 6000N = 64Fs = 2000

1.1. Периодическая последовательность с периодом N/2:

Вывести графики амплитудного и фазового спектра периодической последовательности.

Код программы:

%1C ПЕРИОДИЧЕСКАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ С ПЕРИОДОМN/2

Nb=15;

N=64;

Fs=2000;

T=1/2000;

A1=1.015;

A2=2.03;

f1=250;

f2=6000;

n = 0:(N/2-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ

k = 0:(N/2-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА

w1 = 2*pi*f1/Fs;w2 = 2*pi*f2/Fs;% НОРМИРОВАННЫЕ ЧАСТОТЫ ДИСКРЕТНЫХ ГАРМОНИК

x =A1*cos(w1*n+pi/4)+A2*cos(w2*n+pi/16);% ПЕРИОДИЧЕСКАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ

X =fft(x);% ДПФ ПЕРИОДИЧЕСКОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

MOD = (4/N)*abs(X);% АМПЛИТУДНЫЙ СПЕКТР ПЕРИОДИЧЕСКОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

MOD(1) = (2/N)*abs(X(1));

PHASE =angle(X);% ФАЗОВЫЙ СПЕКТР ПЕРИОДИЧЕСКОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

for i = 1:N/2

if (abs(X(i)) < 1e-4)

      PHASE(i)=0;

end

end

figure('Name','Amplitude Spectrum 1','NumberTitle','off')

subplot(2,1,1), stem(k,MOD,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('1/N|X(k)|')

title(strcat(['Amplitude Spectrum of the Periodic Sequence N/2 = ',num2str(N/2)]))

subplot(2,1,2), stem(k*(2*Fs/N),MOD,'MarkerSize',3,'Linewidth',2),grid

xlabel('f (Hz)'), ylabel('1/N|X(f)|')

title(strcat(['Amplitude Spectrum of the Periodic Sequence N/2 = ',num2str(N/2)]))

figure('Name','Phase Spectrum 1','NumberTitle','off')

subplot(2,1,1), stem(k, PHASE,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('arg{X(k)} (rad)')

title(strcat(['Phase Spectrum of the Periodic Sequence N/2 = ',num2str(N/2)]))

subplot(2,1,2), stem(k*(2*Fs/N),PHASE,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f (Hz)'), ylabel('arg{X(f)} (rad)')

title(strcat(['Phase Spectrum of the Periodic Sequence N/2 = ',num2str(N/2)]))

e1 = 1e-7;% ЗНАЧЕНИЕ ПОРОГА ДЛЯ ПЕРВОГО КРИТЕРИЯ

[MODm,m] =fft_e1(MOD,e1)% ВНЕШНЯЯ ФУНКЦИЯ ДЛЯ ВЫДЕЛЕНИЯ АМПЛИТУД И ЧАСТОТ ГАРМОНИК ПОЛЕЗНОГО СИГНАЛА ПО ПЕРВОМУ КРИТЕРИЮ

Рис. 1 График амплитудного спектра периодической последовательности

Рис. 2 График фазового спектра

Выводы амплитуды и частоты дискретных гармоник:

MODm = 1.9910    1.0150    1.0150 -вектор значений модуля ДПФ, соответствующих сигналу.

m = 0     4    28 -вектор значений частот, соответствующихMODm.

1.2. Конечная последовательность x(n) длины N/2.

Вывести графики модуля и аргумента ДПФ конечной последовательности.

Код программы:

%2С КОНЕЧНАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ ДЛИНЫN/2

MOD_K =abs(fft(x));% МОДУЛЬ ДПФ КОНЕЧНОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

PHASE_K = angle(fft(x));% АРГУМЕНТ ДПФ

for i = 1:N/2

if (abs(X(i)) < 1e-4)

      PHASE_K(i)=0;

end

end

figure('Name','DFT Modulus and Argument 2','NumberTitle','off')

subplot(2,1,1), stem(k,MOD_K,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('|X(k)|')

title('DFT Modulus of the Finite Sequence')

subplot(2,1,2), stem(k, PHASE_K,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('arg{X(k)} (rad)')

title(strcat(['Argument of the Finite Sequence N/2 = ',num2str(N/2)]))

Рис. 3 График модуля и аргумента ДПФ конечной последовательности.

1.3. Конечная последовательность длины N:

Вывести графики конечных последовательностей x(n) и x1(n) и модулей их ДПФ.

Код программы:

%3С КОНЕЧНАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ ДЛИНЫN

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ

k = 0:(N-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА

x0 =A1*cos(w1*n+pi/4)+A2*cos(w2*n+pi/8);% ПЕРИОДИЧЕСКАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ

MOD0 = abs(fft(x0));

for i=1:N/2;

   x1(i) = A1*cos(w1*i);

end

for i=N/2+1:N;

   x1(i) = A2*cos(w2*i);

end

MOD1 = abs(fft(x1));

%ГРАФИК КОНЕЧНЫХ ПОСЛЕДОВАТЕЛЬНОСТЕЙ

figure('Name','Finite Sequence x0','NumberTitle','off')

subplot(2,1,1), stem(n,x0,'MarkerSize',3,'Linewidth',2)

grid, xlabel('n')

ylabel('x(n)'), title(strcat(['Finite Sequence x(n)  N = ',num2str(N)]))

subplot(2,1,2), stem(n,x1,'MarkerSize',3,'Linewidth',2)

grid, xlabel('n')

ylabel('x1(n)'), title(strcat(['Periodic Sequence x1(n)  N = ',num2str(N)]))

%МОДУЛИ ИХ ДПФ

figure('Name','DFT Modulus of Finite Sequence 3','NumberTitle','off')

subplot(2,1,1), stem(k,MOD0,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('|X0(k)|')

title('DFT Modulus x0')

subplot(2,1,2), stem(k,MOD1,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('|X1(k)|')

title('DFTModulusx1')

Рис. 4. График конечных последовательностей

Рис. 5. График модулей их ДПФ.

1.4. Цифровой единичный импульс на интервале n [0; N−1].

Вывести графики цифрового единичного импульса и модуля его ДПФ.

Код программы:

%4С ЦИФРОВОЙ ЕДИНИЧНЫЙ ИМПУЛЬС НА ИНТЕРВАЛЕ [0;N-1]

d1(1)=1;

for i=2:N;

   d1(i) = 0;

end

X4 =fft(d1);

MOD4 =abs(X4);% АМПЛИТУДНЫЙ СПЕКТР ПЕРИОДИЧЕСКОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

figure('Name','Periodic Sequence','NumberTitle','off')

subplot(3,1,1), stem(n,d1,'MarkerSize',3,'Linewidth',2)

grid, xlabel('n')

ylabel('d1(n)'), title(strcat(['Periodic Sequence d1(n)  N = ',num2str(N)]))

subplot(3,1,2), stem(n/Fs,d1,'MarkerSize',3,'Linewidth',2)

grid, xlabel('nT')

ylabel('d1(nT)'), title(strcat(['Periodic Sequence d1(nT)  N = ',num2str(N)]))

d1 =ifft(X4);% ПЕРИОДИЧЕСКАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ, ВЫЧИСЛЕННАЯ С ПОМОЩЬЮ ОДПФ

subplot(3,1,3), stem(n,d1,'MarkerSize',3,'Linewidth',2)

grid, xlabel('n')

ylabel('d1(n)'), title(strcat(['Periodic Sequence d1 = ifft(X)  N = ',num2str(N)]))

figure('Name','DFT Modulus of digital single impulse','NumberTitle','off')

subplot(2,1,1), stem(k,MOD4,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('|X(k)|')

title(strcat(['DFT Modulus of digital single impulse N = ',num2str(N)]))

subplot(2,1,2), stem(k*(Fs/N),MOD4,'MarkerSize',3,'Linewidth',2),grid

xlabel('f (Hz)'), ylabel('|X(f)|')

title(strcat(['DFT Modulus of digital single impulse N = ',num2str(N)]))

Рис. 6. График цифрового единичного импульса

Рис. 7. График модулей их ДПФ.

1.5. Последовательность с однотональной амплитудной модуляцией

Период последовательности 2N.

Вывести графики последовательности и ее амплитудного спектра.

Код программы:

%5C ПОСЛЕДОВАТЕЛЬНОСТЬ С ОДНОТОНАЛЬНОЙ АМПЛИТУДНОЙ МОДУЛЯЦИЕЙ

n = 0:(2*N-1);

k = 0:(2*N-1);

w0 = 2*pi/4;W = w0/4;

fi0 = 0; fiW = 0; m = 0.5;

C = 1;

x5 = C*(cos(w0*n)+m*(0.5*cos((w0-W)*n)+0.5*cos((w0+W)*n)));

X5 = fft(x5);

MOD5 = (2/N)*abs(X5);

MOD5(1) = (1/N)*abs(X5(1));

figure('Name','Periodic Sequence x5','NumberTitle','off')

subplot(2,1,1), stem(n,x5,'MarkerSize',3,'Linewidth',2)

grid, xlabel('n')

ylabel('x5(n)'), title(strcat(['Periodic Sequence x5(n)  N = ',num2str(N)]))

subplot(2,1,2), stem(n/Fs,x5,'MarkerSize',3,'Linewidth',2)

grid, xlabel('nT')

ylabel('x5(nT)'), title(strcat(['Periodic Sequence x(nT)  N = ',num2str(N)]))

figure('Name','Amplitude Spectrum x5','NumberTitle','off')

subplot(2,1,1), stem(k,MOD5,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('1/N|X5(k)|')

title(strcat(['Amplitude Spectrum of the Periodic Sequence N = ',num2str(N)]))

subplot(2,1,2), stem(k*(Fs/N),MOD5,'MarkerSize',3,'Linewidth',2),grid

xlabel('f (Hz)'), ylabel('1/N|X5(f)|')

title(strcat(['Amplitude Spectrum of the Periodic Sequence N = ',num2str(N)]))

Рис. 8. Графикпоследовательности

Рис. 9. Графикее амплитудного спектра.

1.6. Последовательность x(t):

Вывести графики последовательности на интервале t = nT[−500(N−1)T ;500(N−1)T ] с шагом T и модуля ее ДПФ.

Код программы:

%6C ПОСЛЕДОВАТЕЛЬНОСТЬx(t)=sin(pi*n*T)/(pi*n*T)

n = (-500*(N-1)):(500*(N-1));

k = (-500*(N-1)):(500*(N-1));

fd = 2000;

T = 1/fd;

x6 = sinc(pi*T*n);

X6 = fft(x6);

MOD6 = abs(X6);

figure('Name','Periodic Sequence x6','NumberTitle','off')

subplot(2,1,1), stem(n,x6,'MarkerSize',3,'Linewidth',2)

grid, xlabel('n')

ylabel('x6(n)'), title(strcat(['Periodic Sequence x6(n)  N = ',num2str(500*(N-1))]))

subplot(2,1,2), stem(n/fd,x6,'MarkerSize',3,'Linewidth',2)

grid, xlabel('nT')

ylabel('x6(nT)'), title(strcat(['Periodic Sequence x6(nT)  N = ',num2str(500*(N-1))]))

figure('Name','DFT Modulus of sequence x6','NumberTitle','off')

subplot(2,1,1), stem(k,MOD6,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('|X(k)|')

title(strcat(['DFT Modulus of sequence x6 N = ',num2str(500*(N-1))]))

subplot(2,1,2), stem(k*(Fs/N),MOD6,'MarkerSize',3,'Linewidth',2),grid

xlabel('f (Hz)'), ylabel('|X(f)|')

title(strcat(['DFT Modulus of sequence x6 N = ',num2str(500*(N-1))]))

Рис.10. Графикпоследовательности

Рис. 11. Графикее амплитудного спектра.

1.7. Гауссов радиоимпульс:

Задать a = 0,0005 и

Вывести графики последовательности на интервале и модуля ее ДПФ.

Код программы:

%7C Гауссов радиоимпульс

n = (-3*(N-1)):(3*(N-1));

k = (-3*(N-1)):(3*(N-1));

w7 = pi/12;

x7 = exp(-0.0005*n.^2).*cos(w7*n);

X7 = fft(x7);

MOD7 = abs(X7);

figure('Name','Periodic Sequence x7','NumberTitle','off')

subplot(2,1,1), stem(n,x7,'MarkerSize',3,'Linewidth',2)

grid, xlabel('n')

ylabel('x7(n)'), title(strcat(['Periodic Sequence x7(n)  N = ',num2str(3*(N-1))]))

subplot(2,1,2), stem(n/Fs,x7,'MarkerSize',3,'Linewidth',2)

grid, xlabel('nT')

ylabel('x7(nT)'), title(strcat(['Periodic Sequence x7(nT)  N = ',num2str(3*(N-1))]))

figure('Name','DFT Modulus of sequence x7','NumberTitle','off')

subplot(2,1,1), stem(k,MOD7,'MarkerSize',3,'Linewidth',2), grid

xlabel('k'), ylabel('|X7(k)|')

title(strcat(['DFT Modulus of sequence x7 N = ',num2str(3*(N-1))]))

subplot(2,1,2), stem(k*(Fs/N),MOD7,'MarkerSize',3,'Linewidth',2),grid

xlabel('f (Hz)'), ylabel('|X7(f)|')

title(strcat(['DFT Modulus of sequence x7 N = ',num2str(3*(N-1))]))

Рис.12. Графикпоследовательности

Рис. 13. Графикее амплитудного спектра.

Задание 2: Дискретное преобразование Фурье (часть 2)

Исходные данные:

Nb = 15, N = 64, Fs =2000, A1 = 1.015, A2 = 2.03, f1 = 500, f2 = 3000

2.1. Проверку равенства Парсеваля для конечной последовательности с разными длинами N и M.

Код программы:

Nb=15;

N=128;

Fs=2000;

T=1/2000;

A1=1.015;

A2=2.03;

f1=500;

f2=3000;

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ

k = 0:(N-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА

w1 = 2*pi*f1/Fs;w2 = 2*pi*f2/Fs;% НОРМИРОВАННЫЕ ЧАСТОТЫ ДИСКРЕТНЫХ ГАРМОНИК (РАД)

x =A1*cos(w1*n)+A2*cos(w2*n);% ПОСЛЕДОВАТЕЛЬНОСТЬ(ПЕРИОДN)

X =fft(x);% ДПФ ПОСЛЕДОВАТЕЛЬНОСТИ

E1 =sum(x.^2);% ЭНЕРГИЯ ПОСЛЕДОВАТЕЛЬНОСТИ, ВЫЧИСЛЕННАЯ ПО ЕЕ ОТСЧЕТАМ

E2 = (1/N)*sum(abs(X).^2);% ЭНЕРГИЯ ПОСЛЕДОВАТЕЛЬНОСТИ, ВЫЧИСЛЕННАЯ ПО ОТСЧЕТАМ ДПФ

disp(['      E1 = ',num2str(E1),'      E2 = ' num2str(E2)])

Результат проверки:

E1 = 593.4096      E2 = 593.4096

Равенство Парсеваля

2.2. Исследование эффекта растекания спектра

Код программы:

M = 71;% ПЕРИОД ПОСЛЕДОВАТЕЛЬНОСТИM

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ (ПЕРИОДN)

k = 0:(N-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА (ПЕРИОДN)

w1 = 2*pi*f1/Fs;% НОРМИРОВАННАЯ ЧАСТОТА (РАД)

x_N =A1*cos(w1*n);% ПОСЛЕДОВАТЕЛЬНОСТЬ (ПЕРИОДN)

X_N =fft(x_N);% ДПФ ПОСЛЕДОВАТЕЛЬНОСТИ (ПЕРИОДN)

MOD_N = (2/N)*abs(X_N);% АМПЛИТУДНЫЙ СПЕКТР ПОСЛЕДОВАТЕЛЬНОСТИ (ПЕРИОДN)

MOD_N(1) = (1/N)*abs(X_N(1));

n1 = 0:(M-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ (ПЕРИОДM)

k1 = 0:(M-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА (ПЕРИОДM)

x_M =A1*cos(w1*n1);% ПОСЛЕДОВАТЕЛЬНОСТЬ (ПЕРИОДM)

X_M =fft(x_M);% ДПФ ПОСЛЕДОВАТЕЛЬНОСТИ (ПЕРИОДM)

MOD_M = (2/M)*abs(X_M);% АМПЛИТУДНЫЙ СПЕКТР ПОСЛЕДОВАТЕЛЬНОСТИ (ПЕРИОДM)

MOD_M(1) = (1/M)*abs(X_M(1));

P_N =N*f1/Fs;% ЧИСЛО ПЕРИОДОВ ДИСКРЕТНОЙ ГАРМОНИКИ С ЧАСТОТОЙf1 НА ПЕРИОДЕ ПОСЛЕДОВАТЕЛЬНОСТИN

P_M =M*f1/Fs;% ЧИСЛО ПЕРИОДОВ ДИСКРЕТНОЙ ГАРМОНИКИ С ЧАСТОТОЙf1 НА ПЕРИОДЕ ПОСЛЕДОВАТЕЛЬНОСТИM

if P_N-floor(P_N)==0;

  disp(['N = ',num2str(N),'  -->  P_N = ' num2str(P_N)])

disp('%отсутствие эффекта растекания спектра ');

else

   disp(['N = ',num2str(N),'  -->  P_N = ' num2str(P_N)])

   disp('%наличие эффекта растекания спектра ');

win_N  =hamming(N)';% ОКНО ХЭММИНГА — ВЕКТОР-СТОЛБЕЦ ДЛИНЫM

xw_N =x_N.*win_N;% ПОСЛЕДОВАТЕЛЬНОСТЬ, ВЗВЕШЕННАЯ ОКНОМ

XW_N =fft(xw_N);% ДПФ ВЗВЕШЕННОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

MODW_N =(2/N)*abs(XW_N);% АМПЛИТУДНЫЙ СПЕКТР ВЗВЕШЕННОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

MODW_N(1) =(1/N)*abs(XW_N(1));

figure('Name','Reducing Spectrum Leakage with the help of Window Functions','NumberTitle','off')

subplot(2,1,1), stem(k,MOD_N,'MarkerSize',3), grid, xlabel('k')

title(strcat(['Amplitude spectrum without windowing N = ',num2str(N)]))

subplot(2,1,2), stem(k,MODW_N,'MarkerSize',3), grid, xlabel('k')

title(strcat(['Amplitude spectrum with Hamming Window N = ',num2str(N)]))

end

if P_M-floor(P_M)==0;

  disp(['M = ',num2str(M),'  -->  P_M = ' num2str(P_M)])

  disp('%отсутствие эффекта растекания спектра ');

else

   disp(['M = ',num2str(M),'  -->  P_M = ' num2str(P_M)])

   disp('%наличие эффекта растекания спектра ');

win_M  =hamming(M)';% ОКНО ХЭММИНГА — ВЕКТОР-СТОЛБЕЦ ДЛИНЫM

xw_M =x_M.*win_M;% ПОСЛЕДОВАТЕЛЬНОСТЬ, ВЗВЕШЕННАЯ ОКНОМ

XW_M =fft(xw_M);% ДПФ ВЗВЕШЕННОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

MODW_M =(2/M)*abs(XW_M);% АМПЛИТУДНЫЙ СПЕКТР ВЗВЕШЕННОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

MODW_M(1) =(1/M)*abs(XW_M(1));

figure('Name','Reducing Spectrum Leakage with the help of Window Functions','NumberTitle','off')

subplot(2,1,1), stem(k1,MOD_M,'MarkerSize',3), grid, xlabel('k')

title(strcat(['Amplitude spectrum without windowing M = ',num2str(M)]))

subplot(2,1,2), stem(k1,MODW_M,'MarkerSize',3), grid, xlabel('k')

title(strcat(['Amplitude spectrum with Hamming Window M = ',num2str(M)]))

end

Результат исследовании:

N = 128  -->P_N = 32 %отсутствие эффекта растекания спектра

M = 71  -->P_M = 17.75 %наличие эффекта растекания спектра

Рис. 14. Графикамплитудного спектра последовательности до и после применения окна.

P – отношение длительности сигнала к периоду гармонического сигнала

Если Р – целое число то выделение точное, если Р – нецелое число, то наблюдается растекание спектра, потому что вследствие нецелого числа Р в периодическом продолжении с частотой f1 – появляются разрывы (скачки) на границах периода последовательности, из-за которых спектр расширяется.

2.3. Улучшение различения дискретных гармоник с близко расположенными частотами.

Выполнить для конечной последовательностидлины 2N при значениях частоти.

Код программы:

f1_2 = 1.1*f1;% ЧАСТОТЫ ДИСКРЕТНЫХ ГАРМОНИК (Гц)

f2_2 = 1.15*f1;

disp(['      N = ',num2str(N)])

disp(['      f1_2 = ',num2str(f1_2),'       f2_2 = ' num2str(f2_2)])

Delta_N =Fs/N;% РАЗРЕШЕНИЕ ПО ЧАСТОТЕ

Delta_f =abs(f1_2-f2_2);% РАССТОЯНИЕ МЕЖДУ ЧАСТОТАМИ

L =ceil(Fs/(Delta_f-Delta_N));% ВЫБРАННАЯ ДЛИНАL

Delta_L =Fs/L;% ПЕРИОД ДИСКРЕТИЗАЦИИ ПО ЧАСТОТЕ ПРИ ДЛИНЕL

disp(['      Delta_N = ',num2str(Delta_N)])

disp(['      Delta_f = ',num2str(Delta_f)])

disp(['      L = ',num2str(L)])

disp(['      Delta_L = ',num2str(Delta_L)])

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ

w1_2 = 2*pi*f1_2/Fs;w2_2 = 2*pi*f2_2/Fs;% НОРМИРОВАННЫЕ ЧАСТОТЫ

x2 =A1*cos(w1_2*n)+A2*cos(w2_2*n);% КОНЕЧНАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ

X2 =fft(x2);% ДПФ КОНЕЧНОЙ ПОСЛЕДОВАТЕЛЬНОСТИ ДЛИНЫN

MOD2 =abs(X2);% МОДУЛЬ ДПФ

X2_L =fft(x2,L);% ДПФ КОНЕЧНОЙ ПОСЛЕДОВАТЕЛЬНОСТИ, ДОПОЛНЕННОЙ НУЛЯМИ ДО ДЛИНЫL

MOD2_L = abs(X2_L);% МОДУЛЬ ДПФ

k = 0:(N-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА ПРИ ДЛИНЕN

k1 = 0:(L-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА ПРИ ДЛИНЕL

figure('Name','Discrete Harmonic Signal with Close Frequencies','NumberTitle','off')

subplot(2,1,1), stem(k,MOD2), grid, xlabel('k')

title(strcat(['DFT Modulus N = ',num2str(N)]))

subplot(2,1,2), plot(k1,MOD2_L,'r','MarkerSize',3,'Linewidth',2)

grid, holdon, stem(k1,MOD2_L,':'), xlabel('k')

title(strcat(['Spectral Density Modulus L = ',num2str(L)]))

L_2 =ceil(L/2);% ОСНОВНАЯ ПОЛОСА ЧАСТОТL/2

[MODmm]=max(MOD2_L(1:(L_2)));% МАКСИМУМMODm И ИНДЕКСm ВЕКТОРАMOD2_L (ПЕРВЫЙ ПИК)

k_1 = (m-1);f_1 =k_1*Delta_L;% ДИСКРЕТНАЯ НОРМИРОВАННАЯ И АБСОЛЮТНАЯ (Гц) ЧАСТОТЫ ПЕРВОГО ПИКА

K =ceil(L/N);% КОЛИЧЕСТВО ОТСЧЕТОВ НА ПЕРИОДЕ ДИСКРЕТИЗАЦИИFs/N

K1 =m+K;K2 =m+2*K-1;% НИЖНЯЯK1 и ВЕРХНЯЯK2 ГРАНИЦЫ ИНТЕРВАЛА ПРИ ПОИСКЕ ВТОРОГО ПИКА СПРАВА

[MODm1m1]=max(MOD2_L(K1:K2));% МАКСИМУМMODm1 И ИНДЕКСm1 МОДУЛЯ ДПФMOD2_L НА ИНТЕРВАЛЕ [K1K2]

K3 =m-(2*K-1);K4 =m-K;% НИЖНЯЯK3 и ВЕРХНЯЯK4 ГРАНИЦЫ ИНТЕРВАЛА ПРИ ПОИСКЕ ВТОРОГО ПИКА СЛЕВА

[MODm2m2]=max(MOD2_L(K3:K4));% МАКСИМУМMODm2 И ИНДЕКСm2 МОДУЛЯ ДПФMOD2_L НА ИНТЕРВАЛЕ [K3K4]

if (MODm1>MODm2)

k_2 = (K1+m1-1)-1;f_2 =k_2*Delta_L;% ДИСКРЕТНАЯ НОРМИРОВАННАЯ И АБСОЛЮТНАЯ (Гц) ЧАСТОТЫ ВТОРОГО ПИКА, ЕСЛИ ОН СПРАВА ОТ ПЕРВОГО

else

k_2 = (K3+m2-1)-1;f_2 =k_2*Delta_L;% ДИСКРЕТНАЯ НОРМИРОВАННАЯ И АБСОЛЮТНАЯ (Гц) ЧАСТОТЫ ВТОРОГО ПИКА, ЕСЛИ ОН СЛЕВА ОТ ПЕРВОГО

end

disp(['      k_1 = ',num2str(k_1),'      f_1 = ' num2str(f_1)])

disp(['      k_2 = ',num2str(k_2),'      f_2 = ' num2str(f_2)])

Результат вычисления:

N = 128 f1_2 = 550f2_2 = 575

Вывод разрешения по частотеDelta_N, расстояния между частотамиDelta_f, длиныL последовательности и периода дискретизации по частотеDelta_L:

Delta_N = 15.625Delta_f = 25

    L = 214Delta_L = 9.3458

Частоты ближайших пиков:

k_1 = 62f_1 = 579.4393

k_2 = 59      f_2 = 551.4019

Рис. 15. ГрафикN-точечного ДПФ и модуль спектральной плотности, восставленной по L точкам.

Причина погрешности заключается в том, что частотыf12 иf22некратные периоду (Δf = 1/Tд.)

2.4. Вычисление круговой свертки с помощью ДПФ и ОДПФ.

Формула круговой свертки:

ГдеN – длины последовательностейx1 иx2

x(n) – свертка двух последовательностейx1 иx2

Выполнить для последовательностей и

.

Код программы:

x1 = [0 0.25 0.5 0.75 1];% ПЕРВАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ

x2 = [0 0.5 1 0.5 0];% ВТОРАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ

y12 =ifft(fft(x1).*fft(x2));% КРУГОВАЯ СВЕРТКА ПОСЛЕДОВАТЕЛЬНОСТЕЙ

L12 =length(y12);% ПЕРИОД КРУГОВОЙ СВЕРТКИ

figure('Name','Sequences x1, x2, y12','NumberTitle','off')

subplot(3,1,1), stem((0:3*L12-1),...

repmat(x1,1,3),'fill','Linewidth',2,'MarkerSize',3), grid

xlabel('n'), title('Periodic Sequence x1(n)')

subplot(3,1,2), stem((0:3*L12-1), repmat(x2,1,3),'fill','Linewidth',2,'MarkerSize',3), grid

xlabel('n'), title('Periodic Sequence x4(n)')

subplot(3,1,3), stem((0:3*L12-1), repmat(y12,1,3),'fill','Linewidth',2,'MarkerSize',3), grid, xlabel('n')

title('Periodic Sequence y12(n) — Convolution with FFT and IFFT')

Рис. 16. Графикпоследовательностей и их круговой свертки.

2.5. Вычисление реакции ЛДС по формуле свертки с помощью ДПФ.

В качестве воздействия выбрать последовательность с однотональной амплитудной модуляцией:

Период последовательности 2N.

Код программы:

b = [0.8 0.99936 0.64];% КОЭФФИЦИЕНТЫ ЧИСЛИТЕЛЯ ПЕРЕДАТОЧНОЙ ФУНКЦИИ

a = [1 -1.1528 0.73];% КОЭФФИЦИЕНТЫ ЗНАМЕНАТЕЛЯ ПЕРЕДАТОЧНОЙ ФУНКЦИИ

N1 = 25;% ДЛИНА ИМПУЛЬСНОЙ ХАРАКТЕРИСТИКИ

N2 = 128;% ДЛИНА ВОЗДЕЙСТВИЯ

h =impz(b,a,N1)';% ИМПУЛЬСНАЯ ХАРАКТЕРИСТИКА

n = 0:(N2-1);

k = 0:(N2-1);

w0 = 2*pi/4;W = w0/4;

fi0 = 0; fiW = 0; m = 0.5;

C = 1;

x5 = C*(cos(w0*n)+m*(0.5*cos((w0-W)*n)+0.5*cos((w0+W)*n)));% ВОЗДЕЙСТВИЕ

y5_1 =conv(x5,h);% РЕАКЦИЯ, ВЫЧИСЛЕННАЯ С ПОМОЩЬЮ ФУНКЦИИconv

y5_2 =fftfilt(h,x5);% РЕАКЦИЯ, ВЫЧИСЛЕННАЯ С ПОМОЩЬЮ ФУНКЦИИfftfilt

L=N1+N2-1;% ДЛИНА СВЕРТКИ, ВЫЧИСЛЕННОЙ С ПОМОЩЬЮ ФУНКЦИИconv

figure('Name','Impulse Response, Input and Output Signals','NumberTitle','off')

subplot(4,1,1)

stem(0:length(h)-1,h,'Linewidth',2,'MarkerSize',3), grid,

xlabel('n'), title('Impulse Response  h(n)'), xlim([0 L-1])

subplot(4,1,2)

stem(0:length(x5)-1,x5,'Linewidth',2,'MarkerSize',3), grid

xlabel('n'), title('Input Signal x5(n)'), xlim([0 L-1])

subplot(4,1,3)

stem(0:length(y5_1)-1,y5_1,'Linewidth',2,'MarkerSize',3),

grid

xlabel('n'), title('Output Signal y5(n) — Convolution'), xlim([0 L-1])

subplot(4,1,4)

stem(0:length(y5_2)-1,y5_2,'Linewidth',2,'MarkerSize',3), grid

xlabel('n'), title('Output Signal y5(n) — Convolution with FFT and IFFT'), xlim([0 L-1])

При вычислении свертки с помощью ДПФ, количество вычислений сокращается, потому что с помощью алгоритма БПФ, количество интеграций уменьшается вN/log2N раз.

По первому способу (с помощью функцииconv), длина реакции равнаN1 +N2 -1 = 25 + 128 – 1 = 152, соответствует формуле вычисления линейной свертки.

По второму способу (с помощью функцииfftfilt), длина реакции равнаN2 = 128, соответствует формуле вычисления круговой свертки.

При вычислении свертки с помощью ДПФ, длина реакции ограничена (до длины воздействия).

Рис. 17. ГрафикИХ, воздействия и реакции.

2.6. Вычисление реакции ЛДС методом перекрытия с накоплением.

Выполнить с исходными данными п. 2.5 при длине воздействия 4N.

Код программы:

N3 = 4*64;% ДЛИНА ВОЗДЕЙСТВИЯ

x6 =input_1(N3);% ВОЗДЕЙСТВИЕ

y6_1 =fftfilt(h,x6);% РЕАКЦИЯ, ВЫЧИСЛЕННАЯ С ПОМОЩЬЮ ФУНКЦИИfftfilt

y6_2 =fftfilt(h,x6,N1);% РЕАКЦИЯ, ВЫЧИСЛЕННАЯ С ПОМОЩЬЮ ФУНКЦИИfftfilt МЕТОДОМ НАКОПЛЕНИЯ С ПЕРЕКРЫТИЕМ

figure('Name','Impulse Response, Input and Output Signals — Overlap-add method','NumberTitle','off')

subplot(4,1,1)

stem(0:length(h)-1,h,'MarkerSize',3), grid

xlabel('n'), title('Impulse Response h(n)'), xlim([0 N3-1])

subplot(4,1,2), stem(0:length(x6)-1,x6,'MarkerSize',3), grid

xlabel('n'), title('Input Signal x6(n)')

subplot(4,1,3),stem(0:length(y6_1)-1,y6_1,'MarkerSize',3), grid

xlabel('n')

title('Output Signal y6(n) — Convolution with FFT and IFFT')

subplot(4,1,4), stem(0:length(y6_2)-1,y6_2,'MarkerSize',3), grid

xlabel('n')

title('Output Signal y6(n) — Convolution with Overlap-add method')

Рис. 18. ГрафикИХ, воздействия и реакции.

В том случае, когда длина воздействия весьма велика или заранее неизвестна или много больше (чем длина импульсной характеристики) нужно применить метод перекрытия с накоплением или с суммированием.

Задание 3: Спектральный анализ: непараметрические методы

3.1. Проверка оценки СПМ на асимптотическую несмещенность и состоятельность для равномерного белого шума.

Для равномерного белого шума e(n) (идентификатор e) с нулевым средним и единичной дисперсией вычислить:

• истинную СПМ(идентификатор SWN) шума с учетом множителя, задавая длину шума N =100000 (теоретически N → ∞);

• в цикле для каждой длины шума N отNнач =1000 до Nкон =100000 с шагом Δ N = 1000 (вектор N_WN):

периодограмму (идентификатор SWN_estimate);

смещение β (идентификатор beta);

средний квадрат отклонения истинной СПМ от ее оценки (иденти-

фикатор mean_square).

Код программы:

%1C  ПРОВЕРКА ОЦЕНКИ СПМ НА АСИМПТОТИЧЕСКУЮ НЕСМЕЩЕННОСТЬ И СОСТОЯТЕЛЬНОСТЬ

SWN =var(rand(1,100000))/Fs;% ИСТИННАЯ СПМ РАВНОМЕРНОГО БЕЛОГО ШУМА

N_WN = 1000:1000:100000;% ВЕКТОР ДЛИН ШУМА

fori = 1:length(N_WN)% ИНДЕКС ЭЛЕМЕНОВ ВЕКТОРАN_WN

e =rand(1,N_WN(i));% ШУМ ЗАДАННОЙ ДЛИНЫ

SWN_estimate = periodogram(e,[],N_WN(i),Fs,'twosided');% ОЦЕНКА СПМ ШУМА

beta(i) = mean(SWN-SWN_estimate');% СМЕЩЕНИИЕ ОЦЕНКИ СПМ

mean_square (i) =var(SWN_estimate)+beta(i)^2;% СРЕДНИЙ КВАДРАТ ОТКЛОНЕНИЯ ИСТИННОЙ СПМ ОТ ЕЕ ОЦЕНКИ

end

figure('Name','Bias and Mean square Deviation of true PSD from its Estimate','NumberTitle','off')

subplot(2,1,1), plot(N_WN, beta,'LineWidth',2), grid, xlabel('N')

ylabel('beta')

title('Bias')

subplot(2,1,2), plot(N_WN, mean_square,'LineWidth',2), grid, xlabel('N')

ylabel('meanerr')

title('Mean square Deviation of true PSD from its Estimate')

Рис 19. графики зависимости смещения и среднего квадрата отклонения истинной СПМ от длины шума N

Из графики видно, что оценка СПМравномерного белого шума является асимптотически несмещенной (при длине N → ∞ смещение β стремятся к нулю).

3.2. Формирование случайной последовательности длины N ≥1000 (входной параметр function-файла) с требуемой АКФ Ry(m ) :

Ry(m ) = 1+(|m|+1), m = 0, ... , N-1

с последующей фильтрацией последовательности с целью удаления тренда

с помощью БИХ-фильтра ФВЧ Золотарева—Кауэра 5-го порядка с нормиро-

ванной частотой среза WDn = 0.4, максимально допустимым затуханием в ПП

rp = 0.5 и минимально допустимым затуханием в ПЗ rs = 40.

Код программы:

%2С ФОРМИРОВАНИЕ СЛУЧАЙНОЙ ПОСЛЕДОВАТЕЛЬНОСТИ С ТРЕБУЕМОЙ АКФ

N = 1000;% ДЛИНА ШУМА

xw =rand(1,N);% РАВНОМЕРНЫЙ БЕЛЫЙ ШУМ ДЛИНЫN

N02 =var(xw)+(mean(xw)).^2;% КОНСТАНТАN0/2

m = -(N-1):(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ДЛЯ АКФ, ЦЕНТРИРОВАННОЙ ОТНОСИТЕЛЬНОm = 0

Ry_required = 1./(abs(m)+1.);% ТРЕБУЕМАЯ АКФ

L = 2*N-1;% ДЛИНА АКФ

m = 0:L-1;% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ДЛЯ АКФ, ЦЕНТРИРОВАННОЙ ОТНОСИТЕЛЬНОm =N

figure('Name','Required ACF','NumberTitle','off')

stem(m,Ry_required), grid, xlabel('m'), title('Required ACF Ry')

R=100;

Sy = 2*real(fft(Ry_required(N:L),L))-Ry_required (N);% СПМ В L ТОЧКАХ, ВЫЧИСЛЕННАЯ ПО ТРЕБУЕМОЙ АКФ

A =sqrt(real(Sy)./N02);% АЧХ КИХ-ФИЛЬТРА ВL ТОЧКАХ

k = 0:L-1;% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА

F = -k*pi*R/L;% ЛФЧХ КИХ-ФИЛЬТРА ВL ТОЧКАХ

j =sqrt(-1);

H =A.*exp(j*F);% ЧХ КИХ-ФИЛЬТРА ВL ТОЧКАХ

h1 =real(ifft(H));% ИХ КИХ-ФИЛЬТРА НА ПЕРИОДЕL КРУГОВОЙ СВЕРТКИ

y_ACF =fftfilt(h1,xw);% СЛУЧАЙНАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ ДЛИНЫN С ТРЕБУЕМОЙ АКФ НА ВЫХОДЕ КИХ-ФИЛЬТРА

h =h1(1:R+1);% ИХ ДЛИНЫR+1

Ry_estimate =xcorr(y_ACF)./N;% ОЦЕНКА АКФ РЕАКЦИИ КИХ-ФИЛЬТРА, ЦЕНТРИРОВАННАЯ ОТНОСИТЕЛЬНОN

subplot(2,1,1), stem(m,Ry_required), grid, title('Required ACF Ry')

subplot(2,1,2), stem(m,Ry_estimate), grid, xlabel('m')

title(' ACF Estimate of Output Signal - Ry estimate')

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ДЛЯ ВОЗДЕЙСТВИЯ И РЕАКЦИИ КИХ-ФИЛЬТРА

figure('Name','Impulse Response, Input and Output Signals','NumberTitle','off')

subplot(3,1,1), stem(0:R,h), grid, title('Impulse Response h(n)')

subplot(3,1,2), plot(n,xw), grid

title('Input Signal - uniform White Noise')

subplot(3,1,3), plot(n,y_ACF), grid, xlabel('n')

title('Output Signal with Required ACF - y ACF')

[b,a] = ellip(5,0.5,40,0.4,'high');% КОЭФФИЦИЕНТЫ БИХ-ФИЛЬТРА ФВЧ ЗОЛОТАРЕВА-КАУЭРА

y =filter(b,a,y_ACF);% ПОСЛЕДОВАТЕЛЬНОСТЬ НА ВЫХОДЕ БИХ-ФИЛЬТРА ПОСЛЕ УДАЛЕНИЯ ТРЕНДА

figure('Name','Input and Output Signals of Elliptic filter','NumberTitle','off')

subplot(2,1,1), plot(n,y_ACF), grid, title('Input Signal of Elliptic filter — y ACF')

subplot(2,1,2), plot(n,y), grid, xlabel('n'), ylabel('y(n)')

title('Output Signal of Elliptic filter - y')

Рис. 20. График требуемой АКФ

Рис. 21. График случайной последовательности с требуемой АКФ, оценки ее АКФ.

Рис. 22. График случайной последовательности y(n) на выходе БИХ-фильтра.

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

3.3. Расчет немодифицированной периодограммы случайной последовательности y (n) длины N с выводом графика и значений добротности и СКО.

Метод периодограмм заключается в вычислении оценки СПМ  конечной случайной последовательности длиныN:

Код программы:

%3С РАСЧЕТ ПЕРИОДОГРАММЫ

[S,f] =periodogram(y,[],N,Fs,'twosided');% ПЕРИОДОГРАММА СЛУЧАЙНОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

figure('Name','Periodogram of the uniform white Noise','NumberTitle','off')

subplot(2,1,1), plot(f,S,'Linewidth',2), grid

xlabel('f (Hz)'), ylabel('S(f) (Vt/Hz)')

title('Periodogram of the uniform white noise')

subplot(2,1,2), periodogram(y,[],N,Fs,'twosided')

xlabel('f (Hz)'), ylabel('S(f) (dB/Hz)')

title('Periodogram of the uniform white noise')

disp('   СКО периодограммы')

formatlong% ДЛИННЫЙ ФОРМАТ ДЛЯ ВЫВОДА СКО

STD_S =std(S)% СКО ПЕРИОДОГРАММЫ

format% ВОЗРАТ К ИСХОДНОМУ ФОРМАТУ ДЛЯ ВЫВОДА ДОБРОТНОСТЕЙ

disp('   ДОБРОТНОСТЬ периодограммы')

Q_S =mean(S).^2/var(S)% ДОБРОТНОСТЬ ПЕРИОДОГРАММЫ

Результат вычисления:

СКО периодограммы:

STD_S = 1.004566544355291e-04

ДОБРОТНОСТЬ периодограммы:

Q_S = 0.4133

Рис. 23. График периодограммы

3.4. Расчет периодограммы Даньелла случайной последовательности y(n) длины N для произвольного значения K с выводом графика и значений добротности и СКО.

В методе периодограмм Даньелла сглаживание периодограммы достигается за счет ее усреднения по соседним (2K + 1) частотам, симметрично расположенным относительно текущей частоты ωk:

Код программы:

%4C Расчет периодограммы Даньелла

K = [5 10 20];%ВЕКТОРКОЛИЧЕСТВАУСРЕДНЯЕМЫХЧАСТОТ

figure('Name','Daniell Periodograms for the Different Number of Frequency Intervals','NumberTitle','off')

fori = 1:3% ИНДЕКС ЭЛЕМЕНОВ ВЕКТОРАK

S1 = [S(N-K(i)+1:N);S;S(1:K(i))];% ПЕРИОДОГРАММА, ПЕРИОДИЧЕСКИ ПРОДОЛЖЕННАЯ СЛЕВА И СПРАВА НАK ОТСЧЕТОВ

S2 =smooth(S1,K(i));% РЕЗУЛЬТАТ ВЫЧИСЛЕНИЯ СКОЛЬЗЯЩЕГО СРЕДНЕГО (ВЕКТОР ДЛИНЫN+2K(i))

SD(:,i) =S2(K(i)+1:N+K(i));% ПЕРИОДОГРАММА ДАННЬЕЛЛА (ВЕКТОР ДЛИНЫN) ДЛЯ КОЛИЧЕСТВА УСРЕДНЯЕМЫХ ЧАСТОТK(i)

subplot(4,1,i+1), plot(f, SD(:,i),'Linewidth',2), grid

xlabel('f (Hz)'), ylabel(strcat('SD',num2str(i),'(f)'))

title(strcat(['Daniell Periodogram ',num2str(i),'  Frequency Interval 2K+1, K=',num2str(K(i))]))

end

subplot(4,1,1)

plot(f,S,'Linewidth',2), grid, xlabel('f (Hz)'), ylabel('S(f)')

title('Original non-modified periodogram (Vt/Hz)')

disp('   СКО периодограмм ДАНЬЕЛЛА при различном количестве усредняемых частотK')

STD_SD = [K'std(SD)']% КОЛИЧЕСТВО УСРЕДНЯЕМЫХ ЧАСТОТ И СКО ПЕРИОДОГРАММ ДАНЬЕЛЛА

disp('   ДОБРОТНОСТЬ периодограмм ДАНЬЕЛЛА при различном количестве усредняемых частотK')

Q_SD =mean(SD).^2./var(SD);% ДОБРОТНОСТЬ ПЕРИОДОГРАММ ДАНЬЕЛЛА

Q_SD = [K'Q_SD']

Результат вычисления:

СКО при различном количестве усредняемых частот K

STD_SD = 5.0000    0.0001

  10.0000    0.0001

  20.0000    0.0001

  ДОБРОТНОСТЬ при различном количестве усредняемых частот K

Q_SD =    5.0000    0.9938

  10.0000    1.1837

 20.0000    1.3814

Рис. 24. Графики периодограммы и периодограмм Даньелла.

3.5. Расчет периодограммы Бартлетта случайной последовательности y (n) длины N для произвольной длины фрагмента L с выводом графика и значений добротности и СКО.

В методе периодограмм Бартлетта сглаживание периодограммыдостигается за счет разбиения исходной последовательности длины N на неперекрывающиеся фрагменты длины L:

где p — номер фрагмента, P = N/L — их количество

Код программы:

%5С Расчет периодограммы Бартлетта.

L = [10 20 40];% ВЕКТОР ДЛИН ФРАГМЕНТОВ

figure('Name',' Bartlett Periodograms for Different Fragment Lengths','NumberTitle','off')

fori = 1:3% ИНДЕКС ЭЛЕМЕНОВ ВЕКТОРАL

SB(:,i) =pwelch(y,rectwin(L(i)),0,N,Fs,'twosided');% ПЕРИОДОГРАММА БАРТЛЕТТА ДЛЯ ФРАГМЕНТА ДЛИНЫL(i)

subplot(4,1,i+1), plot(f, SB(:,i),'Linewidth',2), grid

xlabel('f (Hz)'), ylabel(strcat('SB',num2str(i),'(f)'))

title(strcat([' Bartlett Periodogram ',num2str(i),'  L =',num2str(L(i))]))

end

subplot(4,1,1)

plot(f,S,'Linewidth',2), grid, xlabel('f (Hz)'), ylabel('S(f)')

title('Original non-modified periodogram (Vt/Hz)')

disp('   СКО периодограмм БАРТЛЕТТА при различной длине фрагментаL')

STD_SB = [L'std(SB)']% ДЛИНЫ ФРАГМЕНТОВ И СКО ПЕРИОДОГРАММ БАРТЛЕТТА

disp('   ДОБРОТНОСТЬ периодограмм БАРТЛЕТТА при различной длине фрагментаL')

Q_SB =mean(SB).^2./var(SB);% ДОБРОТНОСТЬ ПЕРИОДОГРАММ БАРТЛЕТТА

Q_SB = [L'Q_SB']

Результат вычисления:

СКО периодограмм БАРТЛЕТТА при различной длине фрагмента L

STD_SB =10.0000    0.0000

  20.0000    0.0000

  40.0000    0.0001

Добротность периодограмм БАРТЛЕТТА при различной длине фрагмента L

Q_SB = 10.0000    2.4450

20.0000    1.8642

40.0000    1.5982

Рис. 25. Графики периодограммы и периодограмм Бартлетта.

Дисперсия периодограммы Бартлетта обратно пропорциональна количеству фрагментов P, периодограмма Бартлетта более сглаженной (менее осциллирующей), чем периодограмма Даньелла.

3.6. Расчет периодограммы Уэлча случайной последовательности y(n) длины N для произвольной длины фрагмента L и величины перекрытия Q с окном Хэмминга с выводом графика и значений добротности и СКО.

В методе периодограмм Уэлча сглаживание периодограммы достигается за счет разбиения исходной последовательности длины N на перекрывающиеся фрагменты длины L со сглаживающим окном w(n) и величиной перекрытия Q, где Q < L:

где p — номер фрагмента, P — их количество на длине N исходного сигнала, равное:

Код программы:

%6С Расчет периодограммы Уэлча.

num = 1;% ПОРЯДКОВЫЙ НОМЕР ПЕРИОДОГРАММЫ УЭЛЧА

L =ceil([0.1 0.05 0.025].*length(y));% ВЕКТОР ДЛИН ФРАГМЕНТОВ

figure('Name','Welch Periodograms for Different Fragment Lengths and Overlapping','NumberTitle','off')

fori = 1:3% ИНДЕКС ЭЛЕМЕНОВ ВЕКТОРАL

Q =ceil(0.0125*length(y));% ВЕЛИЧИНА ПЕРЕКРЫТИЯ

SW(:,num) =pwelch(y,L(i),Q,N,Fs,'twosided');% ПЕРИОДОГРАММА УЭЛЧА ПРИ ДЛИНЕ ФРАГМЕНТАL(i) И ВЕЛИЧИНЕ ПЕРЕКРЫТИЯQ

subplot(4,1,i+1), plot(f,SW(:,num),'Linewidth',2), grid

xlabel('f (Hz)'), ylabel('S(f)')

title(strcat(['Welch Periodogram Fragment length L = ',num2str(L(i)),'  Overlapping Q = ',num2str(Q)]))

num = num+1;

subplot(4,1,1), plot(f,S,'Linewidth',2), grid, xlabel('f (Hz)')

ylabel('S(f)'), title('Original non-modified periodogram (Vt/Hz)')

end

disp('   СКО периодограммы УЭЛЧА при различной длине фрагментаL и величине перекрытияQ')

LL = [L(1)L(2)L(3)];% ДЛИНЫ ФРАГМЕНТОВ

Q = [QQQ];% ВЕЛИЧИНА ПЕРЕКРЫТИЯ

STD_SW = [LL'Q'std(SW)']% ДЛИНЫ ФРАГМЕНТОВ, ВЕЛИЧИНА ПЕРЕКРЫТИЯ И СКО ПЕРИОДОГРАММ УЭЛЧА

disp('   ДОБРОТНОСТЬ периодограмм УЭЛЧА при различной длине фрагментаL и величине перекрытияQ')

Q_SW =mean(SW).^2./var(SW);% ДОБРОТНОСТЬ ПЕРИОДОГРАММ УЭЛЧА

Q_SW = [LL' Q' Q_SW']

Результат вычисления:

СКО периодограммы УЭЛЧА при различной длине фрагментаL и величине перекрытияQ

STD_SW =100.0000   13.0000    0.0001

  50.0000   13.0000    0.0001

  25.0000   13.0000    0.0000

Добротность периодограмм УЭЛЧА при различной длине фрагментаL и величине перекрытияQ

Q_SW =100.0000   13.0000    1.2478

  50.0000   13.0000    1.4682

25.0000   13.0000    1.7596

Рис. 26. Графики периодограммы и периодограмм Уэлча.

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

3.7. Расчет оценки СПМ по методу Блэкмана—Тьюки случайной последовательности y(n) длины N с заданным окном с выводом графика и значений добротности и СКО.

Оценка СПМ по методу Блэкмана—Тьюки случайной последовательности x n) длины N определяется по формуле:

где-оценка АКФ - четная функция длины L1= 2N1−1 , центрированная относительно m = 0 ; w(m) — весовая функция (окно) длины L1= 2N1−1, центрированная относительно m = 0 .

Код программы:

%7С Расчет оценки СПМ по методу Блэкмана—Тьюки.

N1 =ceil(N/10);% МАКСИМАЛЬНЫЙ СДВИГ ПО ВРЕМЕНИ ДЛЯ ОЦЕНКИ АКФ ОТФИЛЬТРОВАННОГО ШУМА

R1 = (1/N).*xcorr(y,N1);% ОЦЕНКА АКФ ОТФИЛЬТРОВАННОГО ШУМА ДЛИНЫ 2N1+1, ЦЕНТРИРОВАННАЯ ОТНОСИТЕЛЬНОN1+1

R =R1(2:length(R1)-1);% ОЦЕНКА АКФ ОТФИЛЬТРОВАННОГО ШУМА ДЛИНЫ 2N1-1, ЦЕНТРИРОВАННАЯ ОТНОСИТЕЛЬНОN1

L1 =length(R);% ДЛИНА ОЦЕНКИ АКФ

Rw(:,1) =R'.*rectwin(L1);% АКФ, ВЗВЕШЕННАЯ ПРЯМОУГОЛЬНЫМ ОКНОМ

Rw(:,2) =R'.*hamming(L1);% АКФ, ВЗВЕШЕННАЯ ОКНОМ ХЭММИНГА

Rw(:,3) =R'.*bartlett(L1);% АКФ, ВЗВЕШЕННАЯ ОКНОМ БАРТЛЕТТА

name(1).win ='RectangularWindow';% ИМЕНА ОКОН (МАССИВ ЗАПИСЕЙname С ОДНИМ ПОЛЕМwin)

name(2).win ='Hamming Window';

name(3).win ='Bartlett Window';

f = 0:Fs/N:Fs-Fs/N;% ВЕКТОР ЧАСТОТ (Гц)

figure('Name','PSD estimates by the Blackman-Tukey method','NumberTitle','off')

for i = 1:3% НОМЕРА ОКОН

SBT(:,i) = (1/Fs)*(2*real(fft(Rw(N1:L1,i),N))- Rw(N1,i));% ОЦЕНКА СПМ ДЛИНЫ N, ВЫЧИСЛЕННАЯ ПО АКФ, ВЗВЕШЕННОЙ ОКНОМ

subplot(4,1,i+1), plot(f,SBT(:,i),'Linewidth',2), grid

xlabel('f (Hz)'), ylabel(strcat('SBT',num2str(i),'(f)'))

title(['PSD Estimate  -  ',strcat(name(i).win)])

end

subplot(4,1,1)

plot(f,S,'Linewidth',2), grid, xlabel('f (Hz)'), ylabel('S(f)')

title('Original non-modified periodogram (Vt/Hz)')

Рис. 27. Графики периодограммы и оценки СПМ по методу Блэкмана—Тьюки.

3.8. Описание немодифицированной периодограммы последовательности y(n) в виде объекта spectrum с выводом графика.

Код программы:

%8C Описание немодифицированной периодограммы последовательностиy(n)

% в виде объекта spectrum

Hs=spectrum.periodogram;

figure('Name','Spectrum periodogram','NumberTitle','off');

psd(Hs,y,'Fs',Fs);

Рис. 28. Графики немодифицированной периодограммы

После создания объекта spectrum выводится график периодограммы (дБ/Гц) с помощью функции:

psd(Hs,x,'Fs',Fs)

Задание 4: Спектральный анализ: параметрические методы

Исходные данные:

4.1. Моделирование случайной последовательности длины L на основе АР-модели с параметрами, заданными вектором a, и вычисление ее истинной СПМ (Вт/Гц).

АР-модель описывается РУ:

Где - параметры АР-модели, M-1 порядок АР-модели.

АР-модели соответствует БИХ-фильтр полюсного вида ("чисто рекурсивный") с передаточной функцией:

Код программы:

Nb = 15;% НОМЕР ВАРИАНТА

L = 10000;% ДЛИНА ПОСЛЕДОВАТЕЛЬНОСТИ

Fs = 5000;% ЧАСТОТА ДИСКРЕТИЗАЦИИ (Гц)

a = [1 -0.20 0.20 -0.40 0.10 -0.30 0.20 -0.20];% ВЕКТОР ПАРАМЕТРОВ АР-МОДЕЛИ

n = 0:(L-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ

e =randn(1,L);% ВОЗДЕЙСТВИЕ - НОРМАЛЬНЫЙ БЕЛЫЙ ШУМ ДЛИНЫL

y =filter(1,a,e);% МОДЕЛИРУЕМАЯ СЛУЧАЙНАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ ДЛИНЫL

figure('Name','AR-sequence and True PSD','NumberTitle','off')

subplot(2,1,1), plot(n,y), grid, xlabel('n'), ylabel('y(n)')

title('AR-sequence')

f = 0:Fs/(L-1):Fs;% ВЕКТОР ЧАСТОТ (Гц)

H =freqz(1,a,f,Fs);% КОМПЛЕКСНАЯ ЧАСТОТНАЯ ХАРАКТЕРИСТИКА

S = (1/Fs)*abs(H).^2;% ИСТИННАЯ СПМ МОДЕЛИРУЕМОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

subplot(2,1,2),plot(f,S), grid, xlabel('f (Hz)'), ylabel('S(f)')

title('True PSD')

Рис 29. Графики моделируемой последовательности и истинной СПМ.

АР-модель считается наиболее подходящей для оценки СПМ с острыми пиками, но без глубоких впадин. На практике наибольшее распространение получила АР-модель. Это обусловлено ее адекватностью широкому классу реальных сигналов и наименьшими вычислительными затратами.

4.2. Оценка оптимального порядка АР-модели на основе критерия Байеса.

Информационный критерий Байеса:

Код программы:

%2C Оценка оптимального порядка АР-модели на основе критерия Байеса.

x =y;% АНАЛИЗИРУЕМАЯ ПОСЛЕДОВАТЕЛЬНОСТЬ ПОЛАГАЕТСЯ РАВНОЙ МОДЕЛИРУЕМОЙ

pmin = 1;pmax = 3*(length(a)-1);% МИНИМАЛЬНОЕ И МАКСИМАЛЬНОЕ ЗНАЧЕНИЯ ПОРЯДКА

p = [pmin:pmax];% ВЕКТОР ЗНАЧЕНИЙ ПОРЯДКА АР-МОДЕЛИ

fori =pmin:pmax

[aYWD] =aryule(x,p(i));% ПАРАМЕТРЫ АР-МОДЕЛИ И СРЕДНИЙ КВАДРАТ ОШИБКИ ЛИНЕЙНОГО ПРЕДСКАЗАНИЯ

BIC(i) =L*log(D)+p(i)*log(L);% ЗНАЧЕНИЕ КРИТЕРИЯ БАЙЕСА

variance(i) =D;% СРЕДНИЙ КВАДРАТ ОШИБКИ ЛИНЕЙНОГО ПРЕДСКАЗАНИЯ — ОЦЕНКА ДИСПЕРСИИ НОРМАЛЬНОГО БЕЛОГО ШУМА В АР-МОДЕЛИ

end

[BIC_minp_opt] =min(BIC);% МИНИМАЛЬНОЕ ЗНАЧЕНИЕ КРИТЕРИЯ БАЙЕСА И ОПТИМАЛЬНЫЙ ПОРЯДОК МОДЕЛИ

figure('Name','Mean Square of the Linear Prediction Error and Bayesian Information Criterion','NumberTitle','off')

% ГРАФИКИ зависимостей СРЕДНЕГО КВАДРАТА ОШИБКИ ЛИНЕЙНОГО ПРЕДСКАЗАНИЯ и ЗНАЧЕНИЙ КРИТЕРИЯ БАЙЕСА от ПОРЯДКА МОДЕЛИ

subplot(2,1,1), plot(p,variance,'Linewidth',2), grid

xlabel('p'), ylabel('D'), title('Mean Square of the Linear Prediction Error')

subplot(2,1,2), plot(p,BIC,'r','Linewidth',2), grid

xlabel('p'), ylabel('BIC')

title('Bayesian Information Criterion')

Рис. 30.Графики зависимостей среднего квадрата ошибок линейного предсказания и значений критерия Байеса BIC.

По графику определяется оптимальный порядок АР-модели:

4.3. Вычисление оценки СПМ различными методами.

Код программы:

p_opt = 7% ОПТИМАЛЬНЫЙ ПОРЯДОК АР-МОДЕЛИ

[aYWDYW] =aryule(x,p_opt);% ОЦЕНКА ПАРАМЕТРОВ АР-МОДЕЛИ ПО МЕТОДУ ЮЛА-УОЛКЕРА

[aBURGDBURG] =arburg(x,p_opt);% ОЦЕНКА ПАРАМЕТРОВ АР-МОДЕЛИ ПО МЕТОДУ БЕРГА

[aCOVDCOV] =arcov(x,p_opt);% ОЦЕНКА ПАРАМЕТРОВ АР-МОДЕЛИ КОВАРИАЦИОННЫМ МЕТОДОМ

[aMCOVDMCOV] =armcov(x,p_opt);% ОЦЕНКА ПАРАМЕТРОВ АР-МОДЕЛИ МОДИФИЦИРОВАННЫМ КОВАРИАЦИОННЫМ МЕТОДОМ

D =var(e);% ОЦЕНКА ДИСПЕРСИИ НОРМАЛЬНОГО БЕЛОГО ШУМА ДЛИНЫL

disp('% Вывод ИСТИННЫХ ПАРАМЕТРОВ АР-МОДЕЛИ и их ОЦЕНОК')

disp(' ИСТИННЫЕ ПАРАМЕТРЫ');

disp(['  a =    [' num2str(a),']'])

disp(' МЕТОД ЮЛА-УОЛКЕРА');

disp([' aYW =   [' num2str(aYW),']']);

disp(' МЕТОД БЕРГА');

disp([' aBURG = [' num2str(aBURG),']']);

disp('КОВАРИАЦИОННЫЙМЕТОД');

disp([' aCOV =  [' num2str(aCOV),']']);

disp(' МОДИФИЦИРОВАННЫЙ КОВАРИАЦИОННЫЙ МЕТОД');

disp([' aMCOV = [' num2str(aMCOV) ,']']);

disp('% Вывод СРЕДНЕГО КВАДРАТА ОШИБОК ЛИНЕЙНОГО ПРЕДСКАЗАНИЯ ')

disp([' МЕТОД ЮЛА-УОЛКЕРА:DYW = 'num2str(DYW)]);

disp([' МЕТОД БЕРГА:                             DBURG = ' num2str(DBURG)]);

disp([' КОВАРИАЦИОННЫЙ МЕТОД:                    DCOV = ' num2str(DCOV)]);

disp([' МОДИФИЦИРОВАННЫЙ КОВАРИАЦИОННЫЙ МЕТОД:DMCOV = 'num2str(DMCOV)]);

disp('% Вывод ОЦЕНКИ ДИСПЕРСИИ НОРМАЛЬНОГО БЕЛОГО ШУМА')

disp([' ОЦЕНКА ДИСПЕРСИИ НОРМАЛЬНОГО БЕЛОГО ШУМА:D = 'num2str(D)]);

% ПРОВЕРКА УСТОЙЧИВОСТИ БИХ-ФИЛЬТРА

figure('Name','Z-planezero-poleplots','NumberTitle','off')

subplot(2,2,1), zplane(1,aYW), title ('Yule-Walker method'), grid

xlabel('Re'), ylabel('jIm')

subplot(2,2,2), zplane(1,aBURG), title ('Burg method'), grid

xlabel('Re'), ylabel('jIm')

subplot(2,2,3), zplane(1,aCOV), title ('Covariance method'),grid

xlabel('Re'), ylabel('jIm')

subplot(2,2,4), zplane(1,aMCOV),title ('Modified Covariance method')

xlabel('Re'), ylabel('jIm'), grid

% ВЫЧИСЛЕНИЕ ОЦЕНОК СПМ

[SYW,f] = pyulear(x,p_opt,L,Fs,'twosided');% ОЦЕНКА СПМ ПО МЕТОДУ ЮЛА-УОЛКЕРА

[SBURG,f] =  pburg(x,p_opt,L,Fs,'twosided');% ОЦЕНКА СПМ ПО МЕТОДУ БЕРГА

[SCOV,f] =  pcov(x,p_opt,L,Fs,'twosided');% ОЦЕНКА СПМ КОВАРИАЦИОННЫМ МЕТОДОМ

[SMCOV,f] =pmcov(x,p_opt,L,Fs,'twosided');% ОЦЕНКА СПМ МОДИФИЦИРОВАННЫМ КОВАРИАЦИОННЫМ МЕТОДОМ

figure('Name','PSD Estimates','NumberTitle','off')

subplot(4,1,1), plot(f,SYW,'Linewidth',2), grid

xlabel('f (Hz)'), ylabel('SYW(f)'), title(' PSD estimate using Yule-Walker method')

subplot(4,1,2), plot(f,SBURG,'Linewidth',2), grid

xlabel('f (Hz)'), ylabel('SBURG(f)'), title('PSD estimate using Burg method')

subplot(4,1,3), plot(f,SCOV,'Linewidth',2), grid

xlabel('f (Hz)'), ylabel('SCOV(f)'), title('PSD estimate using Covariance method')

subplot(4,1,4), plot(f,SMCOV,'Linewidth',2), grid

xlabel('f (Hz)'), ylabel('SMCOV(f)')

title('PSD estimate using Modified Covariance method')

Значения среднего квадрата ошибок линейного предсказания:

Метод Юла—Уолкера:DYW = 0.99826

Метод Берга:DBURG = 0.99819

Ковариационный метод:DCOV = 0.9985

Модифицированный ковариационный метод:DMCOV = 0.99858

Значение оценки дисперсии нормального белого шума:D = 0.99991

Рис. 31. График карт нулей и полюсов

Рис. 32. График оценок СПМ, вычисленных различными методами.

4. 4. Сравнение оценок СПМ с истинной СПМ.

Для сравнения оценок СПМ с истинной СПМ вывести значения RMSE.

Код программы:

figure('Name','TruePSDandDifferentPSDEstimates','NumberTitle','off')%ГРАФИКИ ОЦЕНОК СПМ и ИСТИННОЙ СПМ

plot(f,S,'Linewidth',2), xlabel('f (Hz)'), ylabel('S(f)'), grid

holdon

plot(f,SYW,'m','Linewidth',2), grid

plot(f,SBURG,'r','Linewidth',2), grid

plot(f,SCOV,'k','Linewidth',2), grid

plot(f,SMCOV,'g','Linewidth',2), grid

legend('True PSD','PSD estimate Yule-Walker method','PSD estimate Burg method','PSD estimate Covariance method','PSD estimate Modified Covariance method', 0);

disp(['% ВЫЧИСЛЕНИЕ ЗНАЧЕНИЙRMSE']);

RMSE_YW =sqrt((1/L).*sum((S'-SYW).^2));%RMSE ДЛЯ ОЦЕНКИ СПМ ПО МЕТОДУ ЮЛА-УОЛКЕРА

RMSE_BURG =sqrt((1/L).*sum((S'-SBURG).^2));%RMSE ДЛЯ ОЦЕНКИ СПМ ПО МЕТОДУ БЕРГА

RMSE_COV =sqrt((1/L).*sum((S'-SCOV).^2));%RMSE ДЛЯ ОЦЕНКИ СПМ КОВАРИАЦИОННЫМ МЕТОДОМ

RMSE_MCOV =sqrt((1/L).*sum((S'-SMCOV).^2));%RMSE ДЛЯ ОЦЕНКИ СПМ МОДИФИЦИРОВАННЫМ КОВАРИАЦИОННЫМ МЕТОДОМ

disp([' МЕТОД ЮЛА-УОЛКЕРА:RMSE = 'num2str(RMSE_YW )]);

disp([' МЕТОД БЕРГА:                          RMSE = ' num2str(RMSE_BURG )]);

disp([' КОВАРИАЦИОННЫЙ МЕТОД:                 RMSE = ' num2str(RMSE_COV)]);

disp([' МОДИФИЦИРОВАННЫЙ КОВАРИАЦИОННЫЙ МЕТОД:RMSE = 'num2str(RMSE_MCOV)]);

Результат вычисления:

Метод Юла—Уолкера: RMSE = 2.2663e-05

Метод Берга: RMSE = 2.2654e-05

Ковариационный метод: RMSE = 2.2731e-05

Модифицированный ковариационный метод: RMSE = 2.2718e-05

Рис. 33. Графики оценок СПМ, вычисленных различными методами, и истинной СПМ.

Метод Юла—Уолкера обычно применяют при анализе длинных последовательностей. Для коротких последовательностей завышенный порядок АР-модели может сопровождаться смещением пиков оценки СПМ и их расщеплением (появлением ложных пиков).

Метод Берга дает удовлетворительные результаты и при анализе коротких последовательностей, однако завышенный порядок АР-модели может также приводить к смещению и расщеплению пиков оценки СПМ.

Ковариационный и модифицированный ковариационный методы обеспечивают более высокую точность при анализе коротких последовательностей, по сравнению с методом Юла—Уолкера с тем же порядком АР-модели. Модифицированный ковариационный метод, как правило, не приводит к расщеплению пиков оценки СПМ и, по сравнению с ковариационным методом, обеспечивает их меньшее смещение.

4.5. Сравнение оценок СПМ с использованием непараметрических и параметрических методов.

Код программы:

[S,f] =periodogram(y,[],L,Fs,'twosided');% ПЕРИОДОГРАММА СЛУЧАЙНОЙ ПОСЛЕДОВАТЕЛЬНОСТИ

figure('Name','PSD Estimates and Periodogram','NumberTitle','off')

subplot(4,1,1), plot(f,S,'Linewidth',2), xlabel('f (Hz)'), ylabel('S(f)'), title('True PSD'), grid

subplot(4,1,2), plot(f,SBURG,'Linewidth',2), grid

xlabel('f (Hz)'), ylabel('SBURG(f)'), title('PSD estimate using Burg method')

subplot(4,1,3), plot(f,S,'Linewidth',2), grid

xlabel('f (Hz)'), ylabel('S(f) (Vt/Hz)'), title('Periodogram of the uniform white noise')

L = 50;% ДЛИНА ФРАГМЕНТОВ

Q =(L/2);% ВЕЛИЧИНА ПЕРЕКРЫТИЯ

[SW,f] =pwelch(y,hamming(length(y)),Q,L,Fs,'twosided');%периодограмма Уэлча  с применением окна Хэмминга

subplot(4,1,4), plot(f,SW,'Linewidth',2), grid

xlabel('f (Hz)'), ylabel('S(f)'), title('Welch Periodogram')

Рис. 34. Графики оценок СПМ и истинной СПМ

4.6. Описание оценок СПМ в виде объектов spectrum и вывести их графики.

Кодпрограммы:

Hs=spectrum.periodogram;

figure('Name','Spectrum periodogram','NumberTitle','off');

psd(Hs,y,'Fs',Fs);

Рис. 34. Графики spectrum периодограммы.

Задание 5: Многоскоростные системы ЦОС

5.1. Моделирование системы однократной интерполяции.

Код программы:

Nb = 15;% НОМЕР БРИГАДЫ

A1 = 1;% АМПЛИТУДЫ ДИСКРЕТНЫХ ГАРМОНИК ВХОДНОГО СИГНАЛА

A2 = 2;

f1 = 25;% ЧАСТОТЫ (Гц) ДИСКРЕТНЫХ ГАРМОНИК ВХОДНОГО СИГНАЛА

f2 = 75;

Ni = 32;% ПЕРИОД ВХОДНОГО СИГНАЛА

Fs_i = 200;% ЧАСТОТА ДИСКРЕТИЗАЦИИ ВХОДНОГО СИГНАЛА

L = 4;% ВЕКТОР КОЭФФИЦИЕНТОВ ИНТЕРПОЛЯЦИИ

disp(['N =Ni = ',num2str(Ni),'Fs =Fs_i = ',num2str(Fs_i)])

disp('%')

disp([' L = [',num2str(L)']'])

disp('%')

N =Ni;% ПЕРИОД ВХОДНОГО СИГНАЛА

Fs =Fs_i;% ЧАСТОТА ДИСКРЕТИЗАЦИИ ВХОДНОГО СИГНАЛА

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ВХОДНОГО СИГНАЛА

x = A1*sin((2*pi*f1/Fs).*n)+A2*sin((2*pi*f2/Fs).*n);% ВХОДНОЙ СИГНАЛ

[y,h]  =interp(x,L);% ВЫХОДНОЙ СИГНАЛ И ИХ КИХ-ФИЛЬТРА (cellarray)

X =fft(x);% ДПФ ВХОДНОГО СИГНАЛА

MODX = (2/N)*abs(X);% АМПЛИТУДНЫЙ СПЕКТР ВХОДНОГО СИГНАЛА

MODX(1)  = (1/N)*abs(X(1));

k = 0:N-1;% ДИСКРЕТНЫЕ НОРМИРОВАННЫЕ ЧАСТОТЫ ВХОДНОГО СИГНАЛА

figure('Name',' Amplitude Spectrums and Magnitude Responses ("new" frequencies) ','NumberTitle','off')

subplot(3,1,1), stem(k*(Fs/N),MODX,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f (Hz)'), title(strcat(['Amplitude Spectrum x(n) N = ',num2str(N)]))

f = 0:Fs/1000:Fs;% СЕТКА ЧАСТОТ ДЛЯ АЧХ КИХ-ФИЛЬТРА В "НОВОЙ" ШКАЛЕ

Y =fft(y);% ДПФ ВЫХОДНОГО СИГНАЛА(cellarray)

Y_D =Y;% ДПФ ВЫХОДНОГО СИГНАЛА(double)

MODY_D = (2/length(Y_D))*abs(Y_D);% АМПЛИТУДНЫЙ СПЕКТР ВЫХОДНОГО СИГНАЛА(double)

MODY_D(1)  = (1/length(Y_D))*abs(Y_D(1));

MODY =MODY_D;% АМПЛИТУДНЫЙ СПЕКТР ВЫХОДНОГО СИГНАЛА(cellarray)

MAG =abs(freqz(h,1,f,Fs));% АЧХ КИХ-ФИЛЬТРА

k_out = 0:length(Y)-1;% ДИСКРЕТНЫЕ НОРМИРОВАННЫЕ ЧАСТОТЫ ВЫХОДНОГО СИГНАЛА

subplot(3,1,2),

stem(k_out*(Fs/length(Y)),MODY,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f ’ (Hz) ("new" frequencies)')

title(strcat(['L = ',num2str(L),' Amplitude Spectrum y(n) and FIR Magnitude Response']))

holdon, plot(f,MAG,'r','Linewidth',2)

f = 0:(Fs*L)/1000:(Fs*L);% СЕТКА ЧАСТОТ ДЛЯ АЧХ КИХ-ФИЛЬТРА В "СТАРОЙ" ШКАЛЕ

MAG =abs(freqz(h,1,f,(Fs*L)));% АЧХ КИХ-ФИЛЬТРА

k_out = 0:length(Y)-1;% ДИСКРЕТНЫЕ НОРМИРОВАННЫЕ ЧАСТОТЫ ВЫХОДНОГО СИГНАЛА

subplot(3,1,3),stem(k_out*(Fs*L/length(Y)),MODY,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f (Hz) ("old" frequencies)'), xlim([0 max(L)*Fs])

title(strcat(['L = ',num2str(L),' Amplitude Spectrum y(n) and FIR Magnitude Response']))

holdon, plot(f,MAG,'r','Linewidth',2)

Рис 35. графики амплитудных спектров и АЧХ КИХ-фильтров.

5.2. Моделирование системы однократной децимации.

Код программы:

Nd = 256;% ПЕРИОД ВХОДНОГО СИГНАЛА

Fs_d = 1600;% ЧАСТОТА ДИСКРЕТИЗАЦИИ ВХОДНОГО СИГНАЛА

M = 4;% ВЕКТОР КОЭФФИЦИЕНТОВ ДЕЦИМАЦИИ

disp(['N =Nd = ',num2str(Nd),'Fs =Fs_d = ',num2str(Fs_d)])

disp([' M = [',num2str(M)']'])

N =Nd;% ПЕРИОД ВХОДНОГО СИГНАЛА

Fs =Fs_d;% ЧАСТОТА ДИСКРЕТИЗАЦИИ ВХОДНОГО СИГНАЛА

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ВХОДНОГО СИГНАЛА

x = A1*sin((2*pi*f1/Fs).*n)+A2*sin((2*pi*f2/Fs).*n);% ВХОДНОЙ СИГНАЛ

y = decimate(x,M,'fir');% ВЫХОДНОЙ СИГНАЛ(cell array)

X =fft(x);% ДПФ ВХОДНОГО СИГНАЛА

MODX = (2/N)*abs(X);% АМПЛИТУДНЫЙ СПЕКТР ВХОДНОГО СИГНАЛА

MODX(1)  = (1/N)*abs(X(1));

k = 0:N-1;% ДИСКРЕТНЫЕ НОРМИРОВАННЫЕ ЧАСТОТЫ ВХОДНОГО СИГНАЛА

figure('Name',' Amplitude Spectrums in Decimation System','NumberTitle','off')

subplot(2,1,1), stem(k*(Fs/N),MODX,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f (Hz)'),title(strcat(['Amplitude Spectrum x(n) N = ',num2str(N)])), xlim([0 Fs])

Y =fft(y);% ДПФ ВЫХОДНОГО СИГНАЛА

Y_D =Y;% ДПФ ВЫХОДНОГО СИГНАЛА(double)

MODY_D = (2/length(Y_D))*abs(Y_D);% АМПЛИТУДНЫЙ СПЕКТР ВЫХОДНОГО СИГНАЛА(double)

MODY_D(1)  = (1/length(Y_D))*abs(Y_D(1));

MODY =MODY_D;% АМПЛИТУДНЫЙ СПЕКТР ВЫХОДНОГО СИГНАЛА(cellarray)

k_out = 0:length(Y)-1;% ДИСКРЕТНЫЕ НОРМИРОВАННЫЕ ЧАСТОТЫ ВЫХОДНОГО СИГНАЛА

subplot(2,1,2),stem(k_out*(Fs/M/length(Y)),MODY,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f (Hz)'), xlim([0 Fs])

title(strcat(['M = ',num2str(M),' Amplitude Spectrum y(n)']))

Рис 36. графики амплитудных спектров входного и выходного сигнала.

5.3. Моделирование систем однократной и многократной интерполяции.

Код программы:

Ni = 32;% ПЕРИОД ВХОДНОГО СИГНАЛА

Fs_i = 200;% ЧАСТОТА ДИСКРЕТИЗАЦИИ ВХОДНОГО СИГНАЛА

L = 8;% ВЕКТОР КОЭФФИЦИЕНТОВ ИНТЕРПОЛЯЦИИ

L1=L/4;

disp([' N = Ni = ',num2str(Ni),' Fs = Fs_i = ',num2str(Fs_i)])

disp('%')

disp([' L = [',num2str(L)']'])

disp('%')

N =Ni;% ПЕРИОД ВХОДНОГО СИГНАЛА

Fs =Fs_i;% ЧАСТОТА ДИСКРЕТИЗАЦИИ ВХОДНОГО СИГНАЛА

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ВХОДНОГО СИГНАЛА

x = A1*sin((2*pi*f1/Fs).*n)+A2*sin((2*pi*f2/Fs).*n);% ВХОДНОЙ СИГНАЛ

[y,h]  =interp(x,L);% ВЫХОДНОЙ СИГНАЛ И ИХ КИХ-ФИЛЬТРА (cellarray)

figure('Name',' Input and Output Signals in Interpolation System','NumberTitle','off')

subplot(3,1,1), stem(n,x), grid, xlabel('n')

title(strcat(['Input Signal x(n) N = ',num2str(N)]))

[y1,h1]  =interp(x,L1);% ВЫХОДНОЙ СИГНАЛ НА ПЕРВОЙ СИСТЕМЕ однократной интерполяции И ИХ КИХ-ФИЛЬТРА (cellarray)

[y2,h2]  =interp(y1,L1);% ВЫХОДНОЙ СИГНАЛ НА ВТОРОЙ СИСТЕМЕ однократной интерполяцииИ ИХ КИХ-ФИЛЬТРА (cellarray)

[y3,h3]  =interp(y2,L1);% ВЫХОДНОЙ СИГНАЛ многократной интерполяции ИХ КИХ-ФИЛЬТРА (cellarray)

subplot(3,1,2), stem(0:length(y)-1,y), grid, xlabel('n')

title(strcat(['L = ',num2str(L),' Output Signal y(n) length(y) =',num2str(length(y))]))

subplot(3,1,3), stem(0:length(y3)-1,y3), grid, xlabel('n')

title(strcat(['L1 = ',num2str(L1),' Output Signal in multiple interpolation system y3(n) length(y3) =',num2str(length(y3))]))

Рис. 37.графики входного и выходных сигналов.

5.4. Моделирование системы однократной передискретизации с полифазной структурой.

Код программы:

N = 128 ;% ПЕРИОД ВХОДНОГО СИГНАЛА

Fs = 800;% ЧАСТОТА ДИСКРЕТИЗАЦИИ ВХОДНОГО СИГНАЛА

L = 3;% КОЭФФИЦИЕНТ ИНТЕРПОЛЯЦИИ

M = 4;% КОЭФФИЦИЕНТ ДЕЦИМАЦИИ

disp(['N =  ',num2str(N),'Fs =  ',num2str(Fs)])

disp('%')

disp([' L/M = ',num2str(L),'/',num2str(M)])

disp('%')

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ВХОДНОГО СИГНАЛА

x = A1*sin((2*pi*f1/Fs).*n)+A2*sin((2*pi*f2/Fs).*n);% ВХОДНОЙ СИГНАЛ

figure('Name',' Input and Output Signals in Resampling System','NumberTitle','off')

subplot(3,1,1), stem(n,x), grid, xlabel('n')

title(strcat(['Input Signal x(n) N = ',num2str(N)]))

Hr =mfilt.firsrc(L,M)% ОПИСАНИЕ ПОЛИФАЗНОЙ СТУКТУРЫ СИСТЕМЫ ОДНОКРАТНОЙ передискретизации

y =filter(Hr,x);% ВЫХОДНОЙ СИГНАЛ

ni_start =ceil(length(y)/2+1);% НАЧАЛО УСТАНОВИВШЕГОСЯ РЕЖИМА

yi =y(ni_start:end);% ВЫХОДНОЙ СИГНАЛ В УСТАНОВИВШЕМСЯ РЕЖИМЕ

ni =ni_start:(ni_start+length(yi)-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ДЛЯ УСТАНОВИВШЕГОСЯ РЕЖИМА

subplot(3,1,2), stem(0:length(y)-1,y), grid, xlabel('n')

title(strcat(['L/M = ',num2str(L),'/',num2str(M),' Output Signal y(n) length(y) = ',num2str(length(y))]))

subplot(3,1,3), stem(ni,yi), grid, xlabel('n')

title(strcat(['Output Signal y(n) Starting point n = ',num2str(ni_start),' length(y) = ',num2str(length(yi))]))

X =fft(x);% ДПФ ВХОДНОГО СИГНАЛА

MODX = (2/N)*abs(X);% АМПЛИТУДНЫЙ СПЕКТР ВХОДНОГО СИГНАЛА

MODX(1)  = (1/N)*abs(X(1));

Yi =fft(yi);% ДПФ ВЫХОДНОГО СИГНАЛА

MODYi = (2/length(Yi))*abs(Yi);% АМПЛИТУДНЫЙ СПЕКТР ВЫХОДНОГО СИГНАЛА

MODYi(1)  = (1/length(Yi))*abs(Yi(1));

k = 0:(N-1);% ДИСКРЕТНЫЕ НОРМИРОВАННЫЕ ЧАСТОТЫ ВХОДНОГО СИГНАЛА

ki = 0:(length(Yi)-1);% ДИСКРЕТНЫЕ НОРМИРОВАННЫЕ ЧАСТОТЫ ВЫХОДНОГО СИГНАЛА

f = 0:Fs*(L/M)/1000:(Fs*(L/M));% СЕТКА ЧАСТОТ ДЛЯ АЧХ КИХ-ФИЛЬТРА

figure('Name',' Amplitude Spectrums and Magnitude Response in Polyphase Structure Interpolation System','NumberTitle','off')

subplot(2,1,1), stem(k*(Fs/N),MODX,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f (Hz)'), xlim([0 L*Fs])

title(strcat(['Amplitude Spectrum x(n) N = ',num2str(N)]))

subplot(2,1,2), stem(ki*(Fs*(L/M)/length(Yi)),MODYi,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f (Hz)'), xlim([0 L/M*Fs])

title(strcat(['L/M = ',num2str(L),'/',num2str(M),' Amplitude Spectrum y(n) and FIR Magnitude Response']))

%L/M=2/5

L = 2;% КОЭФФИЦИЕНТ ИНТЕРПОЛЯЦИИ

M = 5;% КОЭФФИЦИЕНТ ДЕЦИМАЦИИ

disp([' N =  ',num2str(N),' Fs =  ',num2str(Fs)])

disp([' L/M = ',num2str(L),'/',num2str(M)])

n = 0:(N-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ВХОДНОГО СИГНАЛА

x = A1*sin((2*pi*f1/Fs).*n)+A2*sin((2*pi*f2/Fs).*n);% ВХОДНОЙ СИГНАЛ

figure('Name',' Input and Output Signals in Resampling System','NumberTitle','off')

subplot(3,1,1), stem(n,x), grid, xlabel('n')

title(strcat(['Input Signal x(n) N = ',num2str(N)]))

Hr =mfilt.firsrc(L,M)% ОПИСАНИЕ ПОЛИФАЗНОЙ СТУКТУРЫ СИСТЕМЫ ОДНОКРАТНОЙ передискретизации

y =filter(Hr,x);% ВЫХОДНОЙ СИГНАЛ

ni_start =ceil(length(y)/2+1);% НАЧАЛО УСТАНОВИВШЕГОСЯ РЕЖИМА

yi =y(ni_start:end);% ВЫХОДНОЙ СИГНАЛ В УСТАНОВИВШЕМСЯ РЕЖИМЕ

ni =ni_start:(ni_start+length(yi)-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ДЛЯ УСТАНОВИВШЕГОСЯ РЕЖИМА

subplot(3,1,2), stem(0:length(y)-1,y), grid, xlabel('n')

title(strcat(['L/M = ',num2str(L),'/',num2str(M),' Output Signal y(n) length(y) = ',num2str(length(y))]))

subplot(3,1,3), stem(ni,yi), grid, xlabel('n')

title(strcat(['Output Signal y(n) Starting point n = ',num2str(ni_start),' length(y) = ',num2str(length(yi))]))

X =fft(x);% ДПФ ВХОДНОГО СИГНАЛА

MODX = (2/N)*abs(X);% АМПЛИТУДНЫЙ СПЕКТР ВХОДНОГО СИГНАЛА

MODX(1)  = (1/N)*abs(X(1));

Yi =fft(yi);% ДПФ ВЫХОДНОГО СИГНАЛА

MODYi = (2/length(Yi))*abs(Yi);% АМПЛИТУДНЫЙ СПЕКТР ВЫХОДНОГО СИГНАЛА

MODYi(1)  = (1/length(Yi))*abs(Yi(1));

k = 0:(N-1);% ДИСКРЕТНЫЕ НОРМИРОВАННЫЕ ЧАСТОТЫ ВХОДНОГО СИГНАЛА

ki = 0:(length(Yi)-1);% ДИСКРЕТНЫЕ НОРМИРОВАННЫЕ ЧАСТОТЫ ВЫХОДНОГО СИГНАЛА

f = 0:Fs*(L/M)/1000:(Fs*(L/M));% СЕТКА ЧАСТОТ ДЛЯ АЧХ КИХ-ФИЛЬТРА

figure('Name',' Amplitude Spectrums and Magnitude Response in Polyphase Structure Interpolation System','NumberTitle','off')

subplot(2,1,1), stem(k*(Fs/N),MODX,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f (Hz)'), xlim([0 L*Fs])

title(strcat(['Amplitude Spectrum x(n) N = ',num2str(N)]))

subplot(2,1,2), stem(ki*(Fs*(L/M)/length(Yi)),MODYi,'MarkerSize',3,'Linewidth',2)

grid, xlabel('f (Hz)'), xlim([0 L/M*Fs])

title(strcat(['L/M = ',num2str(L),'/',num2str(M),' Amplitude Spectrum y(n) and FIR Magnitude Response']))

Рис.38. График входного и выходных сигналов. (L/M=3/4)

Рис.39. График амплитудных спектров. (L/M=3/4)

Рис.40. График входного и выходных сигналов. (L/M=2/5)

Рис.41. График амплитудных спектров. (L/M=2/5)

При моделировании систем с полифазными структурами в реакции системы будет наблюдаться переходный процесс на интервале дискретного нормированного времени n[0; nнач ], а в ДПФ реакции — эффект растекания спектра. Поэтому сравнивать реакцию с воздействием имеет смысл в установившемся режиме.

Задание 6: Адаптивные фильтры

Исходные данные:

6.1. Оценка импульсной характеристики неизвестной ЛДС.

Код программы:

Nb=15;N=46;L=4000;

n = 0:(L-1);%ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ

n1 = 0:N-1;% ИНТЕРВАЛ ДИСКРЕТНОГО НОРМИРОВАННОГО ВРЕМЕНИ ДЛЯ ОЦЕНОК ИМПУЛЬСНОЙ ХАРАКТЕРИСТИКИ

disp('% ОЦЕНКА ИМПУЛЬСНОЙ ХАРАКТЕРИСТИКИ НЕИЗВЕСТНОЙ ЛДС - БИХ-ФИЛЬТРА ФНЧ (IIR)')

x =rand(1,N);% ВХОДНОЙ СИГНАЛ НЕИЗВЕСТНОЙ ЛДС равномерный белый шум;

n = 0:length(x)-1;

Px =var(x);% СРЕДНИЙ КВАДРАТ ВХОДНОГО СИГНАЛА АФ

mu_max = 2/(N*Px);% МАКСИМАЛЬНЫЙ ШАГ АДАПТАЦИИ

mu = 0.5*mu_max;% ЗАДАННЫЙ ШАГ АДАПТАЦИИ

epsilon = 1e-6;% КОНСТАНТА, ОПРЕДЕЛЯЮЩАЯ МАКСИМАЛЬНЫЙ ШАГ АДАПТАЦИИ

Hnlms =adaptfilt.nlms(N,1,1,epsilon)% СТРУКТУРА АФ С АЛГОРИТМОМNLMS

disp('% Для вывода ГРАФИКОВ СИГНАЛА ОШИБКИ АФ нажмите <ENTER>')

pause

R1 =round(N/2);% ПОРЯДОК БИХ-ФИЛЬТРА ФНЧ

WDn = 0.3;% НОРМИРОВАННАЯ ЧАСТОТА СРЕЗА БИХ-ФИЛЬТРА ФНЧ

[b,a] =butter(R1,WDn,'low');% КОЭФФИЦИЕНТЫ БИХ-ФИЛЬТРА ФНЧ

d =filter(b,a,x);% ВЫХОДНОЙ СИГНАЛ НЕИЗВЕСТНОЙ ЛДС (IIR)

h =impz(b,a,N);% ИСТИННАЯ ИМПУЛЬСНАЯ ХАРАКТЕРИСТИКА НЕИЗВЕСТНОЙ ЛДС (IIR) ДЛИНЫN (ВЕКТОР-СТОЛБЕЦ)

[y_nlms,e_nlms] =filter(Hnlms,x,d);% ВЫХОДНОЙ СИГНАЛ И СИГНАЛ ОШИБКИ АФ С АЛГОРИТМОМNLMS

h_nlms =Hnlms.coefficients;% ПАРАМЕТРЫ АФ С АЛГОРИТМОМLMS

figure('Name','Error signal of AF for NLMS','NumberTitle','off')

plot(n,e_nlms), grid, xlabel('n'), title('Error signal for NLMS')

disp('%')

disp('% Для вывода ГРАФИКОВ ИСТИННОЙ ИМПУЛЬСНОЙ ХАРАКТЕРИСТИКИ (IIR) и ее ОЦЕНОК нажмите <ENTER>')

pause

figure('Name','True Impulse response IIR and its Estimates','NumberTitle','off')

subplot(2,1,1), stem(n1,h), grid, xlabel('n1')

title('True Impulse response IIR - h(n) with length N')

subplot(2,1,2), stem(n1,h_nlms), grid, xlabel('n1')

title('NLMS Impulse response IIR - h nlms(n)')

disp('% вывод СРЕДНЕГО АБСОЛЮТНОГО ОТКЛОНЕНИЯ')

disp('% ОТСЧЕТОВ ИМПУЛЬСНОЙ ХАРАКТЕРИСТИКИ от ее ОЦЕНОК (IIR)')

disp(['         norm1_nlms = ',num2str((1/N)*norm(h'-h_nlms,1))])

Рис.42. График сигнала ошибки АФ.

Рис.43. График истинной ИХ неизвестной ЛДС и ее оценки.

Вывод среднее абсолютное отклонение отсчетов оценки ИХ от их истинных значений: norm1_nlms = 0.057807

6.2. Очистка сигнала от шума.

Код программы:

Fs = 2000;% ЧАСТОТА ДИСКРЕТИЗАЦИИ (Гц)

A1 = 0.03;% АМПЛИТУДЫ ДИСКРЕТНЫХ ГАРМОНИК

A2 = 0.02;

A3 = 0.015;

f1 = 160;% ЧАСТОТЫ (Гц) ДИСКРЕТНЫХ ГАРМОНИК

f2 = 340;

f3 = 480;

wc = 0.5;% НОРМИРОВАННАЯ ЧАСТОТА РАЗРЫВА КИХ-ФИЛЬТРА ФНЧ

r_gauss =randn(1,L);% НОРМАЛЬНЫЙ БЕЛЫЙ ШУМ

x_noise =r_gauss;% ВХОДНОЙ СИГНАЛ АФ - НОРМАЛЬНЫЙ БЕЛЫЙ ШУМ

b =fir1(R1,wc);% КОЭФФИЦИЕНТЫ КИХ-ФИЛЬТРА ФНЧ, ИСКАЖАЮЩЕГО ШУМ

x_noiseNEW =filter(b,1,x_noise);% ИСКАЖЕННЫЙ ШУМ НА ВЫХОДЕ КИХ-ФИЛЬТРА ФНЧ

w1 = 2*pi*f1/Fs;w2 = 2*pi*f2/Fs;w3 = 2*pi*f3/Fs;% НОРМИРОВАННЫЕ ЧАСТОТЫ ДИСКРЕТНЫХ ГАРМОНИК (РАД)

n = 0:(L-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ

s =A1*cos(w1*n)+A2*cos(w2*n)+A3*cos(w3*n);% ПОЛЕЗНЫЙ СИГНАЛ

d =s+x_noiseNEW;% ОБРАЗЦОВЫЙ СИГНАЛ АФ

Hrls =adaptfilt.rls(N);% СТРУКТУРА АФ С АЛГОРИТМОМRLS

[y_rls,e_rls] =filter(Hrls,x_noise,d);% ВЫХОДНОЙ СИГНАЛ И СИГНАЛ ОШИБКИ АФ С АЛГОРИТМОМRLS

disp('% Для вывода ГРАФИКОВ ПОЛЕЗНОГО СИГНАЛА, ЕГО АДДИТИВНОЙ СМЕСИ С ШУМОМ')

disp('% и ОЦЕНОК ПОЛЕЗНОГО СИГНАЛАRLS нажмите <ENTER>')

pause

figure('Name','Harmonic Signal, Mixture of Signal and Noise, Estimates of Harmonic Signal (RLS)','NumberTitle','off')

subplot(3,1,1), plot(n,s), grid, title('Harmonic Signal - s(n)')

subplot(3,1,2), plot(n,d), grid, title('Mixture of Harmonic Signal and Noise - d(n)')

subplot(3,1,3), plot(n,e_rls), grid, xlabel('n'), title('Estimate Harmonic Signal (RLS)')

disp('% Для вывода ГРАФИКОВ АМПЛИТУДНЫХ СПЕКТРОВ ПОЛЕЗНОГО СИГНАЛА, ЕГО АДДИТИВНОЙ СМЕСИ С ШУМОМ')

disp('% и ОЦЕНОК ПОЛЕЗНОГО СИГНАЛАRLS')

disp('% в УСТАНОВИВШЕМСЯ РЕЖИМЕ нажмите <ENTER>')

pause

k = 0:(L-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА

f = (k/L)*Fs;% ЧАСТОТА (Гц)

MOD_s = (2/L)*abs(fft(s));% АМПЛИТУДНЫЙ СПЕКТР ПОЛЕЗНОГО СИГНАЛА

MOD_s(1) = (1/L)*abs(fft(s(1)));

MOD_d = (2/L)*abs(fft(d));% АМПЛИТУДНЫЙ СПЕКТР СМЕСИ ПОЛЕЗНОГО СИГНАЛА С ШУМОМ

MOD_d(1) = (1/L)*abs(fft(d(1)));

n_start =round(0.05*L);% НАЧАЛО УСТАНОВИВШЕГОСЯ РЕЖИМА

Ls = (L-n_start)+1;% ДЛИНА ОЦЕНОК СИГНАЛА И ШУМА В УСТАНОВИВШЕМСЯ РЕЖИМЕ

e1_rls =e_rls(n_start:end);% ОЦЕНКА ПОЛЕЗНОГО СИГНАЛА В УСТАНОВИВШЕМСЯ РЕЖИМЕ (RLS)

ks = 0:(Ls-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА ДЛЯ УСТАНОВИВШЕГОСЯ РЕЖИМА

fs = (ks/Ls)*Fs;% ЧАСТОТА (Гц)

MOD_rls = (2/Ls)*abs(fft(e1_rls));% АМПЛИТУДНЫЙ СПЕКТР ОЦЕНКИ ПОЛЕЗНОГО СИГНАЛА В УСТАНОВИВШЕМСЯ РЕЖИМЕ (RLS)

MOD_rls(1) = (1/Ls)*abs(fft(e1_rls(1)));

figure('Name','Amplitude Spectrums of Harmonic Signal, Mixture of Signal and Noise, and Estimates of Harmonic Signal (RLS)','NumberTitle','off')

subplot(3,1,1), stem(f,MOD_s), grid, title('Amplitude Spectrum of Harmonic Signal')

subplot(3,1,2), stem(f,MOD_d), grid, title('Amplitude Spectrum of Harmonic Signal and Noise - d(n)')

subplot(3,1,3), stem(fs,MOD_rls), grid, xlabel('f'), title('Amplitude Spectrum of Estimate Harmonic (RLS)')

ns = 0:(Ls-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ДЛЯ УСТАНОВИВШЕГОСЯ РЕЖИМА

s1 =s(n_start:end);% ПОЛЕЗНЫЙ СИГНАЛ НА ИНТЕРВАЛЕ ВРЕМЕНИ ДЛЯ ЕГО ОЦЕНОК

RMSE_rls =sqrt((1/Ls).*sum((s1-e1_rls).^2));%RMSE ДЛЯ ОЦЕНКИ ПОЛЕЗНОГО СИГНАЛА (RLS)

disp('% Для вывода ЗНАЧЕНИЙRMSE нажмите <ENTER>')

pause

disp(['         RMSE_rls =     ',num2str(RMSE_rls)])

Рис.44. График полезного сигнала, его аддитивной смеси.

Рис.45. График амплитудных спектров данных сигналов.

Вывод значения RMSE:RMSE_rls =     0.0047357

6.3. Выравнивание частотной характеристики неизвестной ЛДС.

Код программы:

A1 = 0.03;A2 = 0.02;A3 = 0.015;% АМПЛИТУДЫ ДИСКРЕТНЫХ ГАРМОНИК

w1 =pi/4;% НОРМИРОВАННЫЕ ДИСКРЕТНЫЕ ЧАСТОТЫ (РАД)

w2 = 3*pi/4;

w3 = 3.5*pi/4;

n = 0:(L-1);% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ

s =A1*cos(w1*n)+A2*cos(w2*n)+A3*cos(w3*n)+4*A1*r_gauss;% ВХОДНОЙ СИГНАЛ НЕИЗВЕСТНОЙ ЛДС - ОБРАЗЦОВЫЙ СИГНАЛ АФ

k = 0:(L-1);% ДИСКРЕТНАЯ НОРМИРОВАННАЯ ЧАСТОТА

f =k*(Fs/L);% ЧАСТОТА (Гц)

S =fft(s);% ДПФ ВХОДНОГО СИГНАЛА НЕИЗВЕСТНОЙ ЛДС

MODS = (2/L)*abs(S);% АМПЛИТУДНЫЙ СПЕКТР ВХОДНОГО СИГНАЛА НЕИЗВЕСТНОЙ ЛДС

MODS(1) = (1/L)*abs(S(1));

disp('% Для вывода ГРАФИКОВ ВХОДНОГО СИГНАЛА НЕИЗВЕСТНОЙ ЛДС и его АМПЛИТУДНОГО СПЕКТРА   нажмите <ENTER>')

pause

figure('Name','Input Signal of Unknown LDS and its Amplitude Spectrum','NumberTitle','off')

subplot(2,1,1), plot(n,s), grid, xlabel('n'), title('Input Signal of Unknown LDS')

subplot(2,1,2), stem(f,MODS), grid, xlabel('f'), title('Amplitude Spectrum')

d =s;% ОБРАЗЦОВЫЙ СИГНАЛ АФ - ВХОДНОЙ СИГНАЛ НЕИЗВЕСТНОЙ ЛДС

D =floor(N/2);% ВЕЛИЧИНА ЗАДЕРЖКИ ОБРАЗЦОГО СИГНАЛА

D_zeros =zeros(1,round(N/2));% ФОРМИРОВАНИЕ НАЧАЛЬНЫХ НУЛЕВЫХ ЗНАЧЕНИЙ ЗАДЕРЖАННОГО ОБРАЗЦОВОГО СИГНАЛА

d_delay = [D_zerosd(1:(length(d)-length(D_zeros)))];% ЗАДЕРЖАННЫЙ ОБРАЗЦОВЫЙ СИГНАЛ

R2 =round(N/7);% ПОРЯДОК КИХ-ФИЛЬТРА ФНЧ

wc = 0.5;% НОРМИРОВАННАЯ ЧАСТОТА РАЗРЫВА КИХ-ФИЛЬТРА ФНЧ

b =fir1(R2,wc);% ИМПУЛЬСНАЯ ХАРАКТЕРИСТИКА НЕИЗВЕСТНОЙ ЛДС

x =filter(b,1,d);% ВЫХОДНОЙ СИГНАЛ НЕИЗВЕСТНОЙ ЛДС - ВХОДНОЙ СИГНАЛ АФ

X =fft(x);% ДПФ ВЫХОДНОГО СИГНАЛА НЕИЗВЕСТНОЙ ЛДС

MODX = (2/L)*abs(X);% АМПЛИТУДНЫЙ СПЕКТР ВЫХОДНОГО СИГНАЛА НЕИЗВЕСТНОЙ ЛДС

MODX(1) = (1/L)*abs(X(1));

disp('% Для вывода ГРАФИКОВ ВЫХОДНОГО СИГНАЛА НЕИЗВЕСТНОЙ ЛДС и его АМПЛИТУДНОГО СПЕКТРА нажмите <ENTER>')

pause

figure('Name','Output Signal of Unknown LDS and its Amplitude Spectrum','NumberTitle','off')

subplot(2,1,1), plot(n,x), grid, xlabel('n'), title('Output Signal of Unknown LDS')

subplot(2,1,2), stem(f,MODX), grid, xlabel('f'), title('Amplitude Spectrum')

Hrls = adaptfilt.rls(N);% СТРУКТУРА АФ С АЛГОРИТМОМ RLS

[y_rls,e_rls] = filter(Hrls,x,d_delay);% ВЫХОДНОЙ СИГНАЛ И СИГНАЛ ОШИБКИ АФ С АЛГОРИТМОМ RLS

h_rls =Hrls.coefficients;% ПАРАМЕТРЫ АФ (ИМПУЛЬСНАЯ ХАРАКТЕРИСТИКА ЕГО КИХ-ФИЛЬТРА)

h_conv =conv(h_rls,b);% ИМПУЛЬСНАЯ ХАРАКТЕРИСТИКА КАСКАДНОГО СОЕДИНЕНИЯ НЕИЗВЕСТНОЙ ЛДС И АФ

n2 = 0:length(b)-1;% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ДЛЯ ИМПУЛЬСНОЙ ХАРАКТЕРИСТИКИ НЕИЗВЕСТНОЙ ЛДС (КИХ-ФИЛЬТРА ФНЧ)

n3 = 0:length(h_conv)-1;% ДИСКРЕТНОЕ НОРМИРОВАННОЕ ВРЕМЯ ДЛЯ СВЕРТКИ ИМПУЛЬСНОЙ ХАРАКТЕРИСТИКИ НЕИЗВЕСТНОЙ ЛДС И АФ

disp('% Для вывода ГРАФИКОВ ИМПУЛЬСНЫХ ХАРАКТЕРИСТИК НЕИЗВЕСТНОЙ ЛДС, КИХ-ФИЛЬТРА В СОСТАВЕ АФ')

disp('% и ИХ КАСКАДНОГО СОЕДИНЕНИЯ  нажмите <ENTER>')

pause

figure('Name','Impulse Response of AF, Unknown LDS, and Cascade Connection','NumberTitle','off')

subplot(3,1,1), stem(n2,b), grid, title('Impulse Response of Unknown LDS')

subplot(3,1,2), stem(n1,h_rls), grid, title('Impulse Response of AF')

subplot(3,1,3), stem(n3,h_conv), grid, xlabel('n'), title('Impulse Response of Cascade Connection')

fa = 0:Fs/500:(Fs/2);% ВЕКТОР ЧАСТОТ ДЛЯ АЧХ

MAG_US =abs(freqz(b,1,fa,Fs));% АЧХ НЕИЗВЕСТНОЙ ЛДС

MAG_AF =abs(freqz(h_rls,1,fa,Fs));% АЧХ КИХ-ФИЛЬТРА В СОСТАВЕ АФ

MAG =MAG_US.*MAG_AF;% АЧХ КАСКАДНОГО СОЕДИНЕНИЯ НЕИЗВЕСТНОЙ ЛДС И АФ

disp('% Для вывода ГРАФИКОВ АЧХ НЕИЗВЕСТНОЙ ЛДС, КИХ-ФИЛЬТРА В СОСТАВЕ АФ,')

disp('% и ИХ КАСКАДНОГО СОЕДИНЕНИЯ нажмите <ENTER>')

pause

figure('Name','Magnitude Response of Unknown LDS, AF, and Cascade Connection','NumberTitle','off')

subplot(3,1,1), plot(fa,MAG_US), grid, title('Magnitude Response of Unknown LDS')

subplot(3,1,2), plot(fa,MAG_AF), grid, title('Magnitude Response of AF')

subplot(3,1,3), plot(fa,MAG), grid, xlabel('f'), title('Magnitude Response of Cascade Connection')

PH_US = angle(freqz(b,1,fa,Fs));% ФЧХ НЕИЗВЕСТНОЙ ЛДС

PH_AF = angle(freqz(h_rls,1,fa,Fs));% ФЧХ КИХ-ФИЛЬТРА В СОСТАВЕ АФ

PH =PH_US+PH_AF;% ФЧХ КАСКАДНОГО СОЕДИНЕНИЯ НЕИЗВЕСТНОЙ ЛДС И АФ

disp('% Для вывода ГРАФИКОВ ФЧХ НЕИЗВЕСТНОЙ ЛДС, КИХ-ФИЛЬТРА В СОСТАВЕ АФ')

disp('% и ИХ КАСКАДНОГО СОЕДИНЕНИЯ нажмите <ENTER>')

pause

figure('Name','Phase Response of Unknown LDS, AF, and Cascade Connection','NumberTitle','off')

subplot(3,1,1), plot(fa,PH_US), grid, title('Phase Response of Unknown LDS')

subplot(3,1,2), plot(fa,PH_AF), grid, title('Phase Response of AF')

subplot(3,1,3), plot(fa,PH), grid, xlabel('f'), title('Phase Response of Cascade Connection')

Y =fft(y_rls);% ДПФ ВЫХОДНОГО СИГНАЛА АФ

MODY = (2/L)*abs(Y);% АМПЛИТУДНЫЙ СПЕКТР ВЫХОДНОГО СИГНАЛА АФ

MODY(1) = (1/L)*abs(Y(1));

disp('% Для вывода ГРАФИКОВ ВЫХОДНОГО СИГНАЛА АФ и его АМПЛИТУДНОГО СПЕКТРА нажмите <ENTER>')

pause

figure('Name','Output Signal and Amplitude Spectrum of AF','NumberTitle','off')

subplot(2,1,1), plot(n,y_rls), grid, xlabel('n'), title('Output Signal of AF')

subplot(2,1,2), stem(f,MODY), grid, xlabel('f'), title('Amplitude Spectrum')

Рис.46. График входного сигнала неизвестной ЛДС и его амплитудного спектра

Рис.47. График выходного сигнала неизвестной ЛДС и его амплитудного спектра.

Рис.48. График импульсных характеристик.

Рис.49. График АЧХ неизвестной ЛДС, КИХ-фильтра в составе АФ и их каскадного соединения.

Рис.50. Графики ФЧХ.

Рис.51. График выходного сигнала АФ и его амплитудного спектра.

6.4. Вычисление оценок параметров АР-модели и оценок параметров линейного предсказания.

Код программы:

a = [1 -0.20 0.20 -0.40 0.10 -0.30 0.20 -0.20];% ВЕКТОР ПАРАМЕТРОВ АР-МОДЕЛИ

e_AR =r_gauss;% ВХОДНОЙ СИГНАЛ АР-МОДЕЛИ

y_AR =filter(1,a,e_AR);% ВЫХОДНОЙ СИГНАЛ АР-МОДЕЛИ

x =y_AR;% АНАЛИЗИРУЕМЫЙ СИГНАЛ ПОЛАГАЕТСЯ РАВНЫМ МОДЕЛИРУЕМОМУ СИГНАЛУ

p =length(a)-1;% ПОРЯДОК АР-МОДЕЛИ

[aYWDYW] =aryule(x,p);% ОЦЕНКА ПАРАМЕТРОВ АР-МОДЕЛИ ПО МЕТОДУ ЮЛА-УОЛКЕРА

d =x;% ОБРАЗЦОВЫЙ СИГНАЛ АФ

x_delay = [0x(1:(L-1))];% ВХОДНОЙ СИГНАЛ АФ - ЗАДЕРЖАННЫЙ ОБРАЗЦОВЫЙ СИГНАЛ

[y_rls,e_rls] =filter(Hrls,x_delay,d);% ВЫХОДНОЙ СИГНАЛ И СИГНАЛ ОШИБКИ АФ С АЛГОРИТМОМRLS

h_rls = -Hrls.coefficients;% ВЕКТОР ОЦЕНОК ПАРАМЕТРОВ ЛИНЕЙНОГО ПРЕДСКАЗАНИЯ

disp('% Для вывода ГРАФИКОВ ВЕКТОРА ЗАДАННЫХ ПАРАМЕТРОВ АР-МОДЕЛИ, их ОЦЕНОК')

disp('% и ВЕКТОРА ОЦЕНОК ПАРАМЕТРОВ ЛИНЕЙНОГО ПРЕДСКАЗАНИЯ нажмите <ENTER>')

pause

figure('Name','True AR parameters, Estimated AR parameters, and Estimated Linear Prediction parameters','NumberTitle','off')

YMAX = max([max(abs(a)) max(abs(aYW)) max(abs(h_rls))]);% МАКСИМАЛЬНЫЙ ПО МОДУЛЮ ЭЛЕМЕНТ ВЕКТОРОВ

subplot(3,1,1), stem(a(2:end)), grid, xlim([0 N-1]), ylim([-YMAX YMAX])

title('True parameters - a')

subplot(3,1,2), stem(aYW(2:end)), grid, xlim([0 N-1]), ylim([-YMAX YMAX])

title('Estimated parameters AR - aYW')

subplot(3,1,3), stem(h_rls), grid, xlim([0 N-1]), ylim([-YMAX YMAX])

title('Estimated Linear Prediction parameters - h rls')

disp('% Для вывода ЗНАЧЕНИЙMAE нажмите <ENTER>')

pause

a_N = [a(2:end)zeros(1,(N-length(a(2:end))))];% ВЕКТОР ИСТИННЫХ ПАРАМЕТРОВ АР-МОДЕЛИ, ДОПОЛНЕННЫЙ НУЛЯМИ ДО ДЛИНЫN

aYW_N = [aYW(2:end)zeros(1,(N-length(aYW(2:end))))];% ВЕКТОР ОЦЕНОК ПАРАМЕТРОВ АР-МОДЕЛИ, ДОПОЛНЕННЫЙ НУЛЯМИ ДО ДЛИНЫN

MAE_AR = (1/N).*sum(abs(a_N-aYW_N));%MAE ОЦЕНОК ПАРАМЕТРОВ АР-МОДЕЛИ

MAE_LP = (1/N).*sum(abs(h_rls-a_N));%MAE ОЦЕНОК ПАРАМЕТРОВ ЛИНЕЙНОГО ПРЕДСКАЗАНИЯ

disp('%')

disp(['         MAE_AR =     ',num2str(MAE_AR)])

disp(['         MAE_LP =     ',num2str(MAE_LP)])

Вывод значения MAE:

MAE_AR =     0.002014

MAE_LP =     0.012663

Рис. 52. Графики вектора заданных параметров АР-модели (без учета a0 =1) и вектора оценок параметров АР-модели и вектора оценок параметров линейного предсказания.

Заключение

В настоящее время к общепризнанным универсальным мировым стандартам в области компьютерных технологий относится программная среда (система) MATLAB, предназначенная для моделирования в самых разных областях науки и техники, в первую очередь, ЦОС.

Широкое распространение MATLAB обусловлено следующими основными достоинствами этой системы:

алгоритмическим языком "сверхвысокого" уровня за счет матричной обработки данных;

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

огромным разнообразием графических средств;

удобными средствами создания и отладки программ;

широким набором программных средств общего (ядро MATLAB) и специального (пакеты расширения Toolbox) назначения;

наличием разнообразных средств GUI (Graphical User Interface — графический интерфейс пользователя) без использования алгоритмического языка в явном виде;

широким набором средств Simulink общего (ядро Simulink) и специального назначения для блочного моделирования динамических систем.

Список использованной литературы
  1. Цифровая обработка сигналов. Моделирование вMatlab. А. И.Солонина, С.М. Арбузов СПб.: БХВ-Петербург, 2008
  2. Цифровая обработка сигналов и MATLAB: учеб. пособие / А. И.Солонина,Д. М. Клионский, Т. В. Меркучева, С. Н. Перов. — СПб.: БХВ-Петербург,2013. — 512 с.: ил.— (Учебная литература для вузов)




Возможно эти работы будут Вам интересны.

1. ЦИФРОВАЯ СИСТЕМА КОММУТАЦИИ 5ESS

2. ЦИФРОВАЯ СИСТЕМА КОММУТАЦИИ Alcatel 1000 E-10

3. Цифровая система коммутации типа EWSD

4. Генераторы электрических сигналов

5. Основы анализа сигналов

6. СОГЛАСОВАННАЯ ФИЛЬТРАЦИЯ ДЕТЕРМИНИРОВАННЫХ СИГНАЛОВ

7. Спектральный анализ случайных сигналов

8. Энергетические характеристики вещественных сигналов

9. ВОЗДЕЙСТВИЕ АМ СИГНАЛОВ НА РЕЗОНАНСНУЮ ЦЕПЬ

10. Математические модели радиотехнических сигналов