Моделювання роботи ЛДС з використанням різницевого рівняння

Національний технічний  університет України

«Київський політехнічний  інститут»

ММІФ

Кафедра МКТМ

 

 

 

 

 

Лабораторна робота № 3

з дисципліни «Теорія сигналів»

«Моделювання роботи ЛДС з використанням різницевого рівняння»

 

 

 

 

 

 

 

 

 

        Виконала:

         студентка гр. ІМ-11

         Ковальчук О.

 

         Перевірила:

         Порєва Г.С.

 

 

 

 

 

 

Київ 2013

Мета роботи: дослідити роботу лінійних систем обробки дискретних сигналів; набути навичок моделювання лінійних стаціонарних дискретних систем в MatLAB.

 

 

Основні теоретичні відомості

Математичне моделювання обробки сигналів лінійною дискретною

системою (ЛДС) включає в себе:

– розрахунок характеристик ЛДС в часовій області, z-області, та частотній області;

– розрахунок реакції ЛДС на вхідний сигнал по різницевому рівнянню;

– аналіз вхідних сигналів та реакцій ЛДС в часовій та частотних областях.

В MatLAB математичною моделлю ЛДС називають співвідношення між входом та виходом у вигляді різницевого рівняння або системи рівнянь, які дозволяють розрахувати реакцію системи на заданий вхідний вплив.

В часовій області основною характеристикою ЛДС є імпульсна характеристика а моделювання роботи ЛДС (розрахунок реакції) виконується

з використанням одного з таких співвідношень між входом та виходом:

1. Різницевого рівняння

2. Формули згортки

3. Системи рівнянь змінних стану.

В частотній області основною характеристикою ЛДС є комплексна частотна характеристика, а також її модуль (амплітудно-частотна характеристика, АЧХ) та аргумент (фазо-частотна характеристика, ФЧХ):

 

Порядок роботи

Моделювання роботи ЛДС з використанням різницевого рівняння.

 

1. . На основі власної дати народження записати два різницевих рівняння:

 (1)

 

clc;

close all;

clear all;

D1=2;

D2=7;

M1=0;

M2=9;

P1=1;

P2=9;

P3=9;

P4=3;

a1=[1]

b1=[M2/10,-(P3-D1)/20,(M1-M2)/20,-P4/30,D2/20,-M2/20]

a2=[2,-(D1+D2)/30,(P4-D2)/40,-D1/50,(M2-D2)/50]

b2=[M2/10,-(P3-D1)/20,(M1-M2)/20,-P4/30,D2/20,-M2/20]

a1 =

     1

b1 =

    0.9000   -0.3500   -0.4500   -0.1000    0.3500   -0.4500

a2 =

    2.0000   -0.3000   -0.1000   -0.0400    0.0400

b2 =

    0.9000   -0.3500   -0.4500   -0.1000    0.3500   -0.4500

2. Сформувати відліки синусоїдального сигналу частоти 10 Гц тривалістю 2 сек. амплітуди 1 В, дискретизованого з частотою 256 Гц. Розрахувати реакцію систем (1) та (2) на отриманий сигнал (функція filter) для двох випадків:

2.1. нульові початкові умови

2.2. випадкові початкові умови  (скористатися функцією rand)

Побудувати графіки вхідного та вихідного сигналів в одному вікні, позначивши точки графіку, що відповідають відлікам, та огинаючі графіків. Побудувати в окремому вікні та порівняти перші 0.1 сек. вхідного та вихідного сигналу. Зробити висновки.

fs=256;

t=0:1/fs:2;

t1=0:1/fs:0.1;

fun=sin(10*t);

y1=filter(b1,a1,fun);

plot(t,y1),hold on;

plot(t,fun);

title('Графік sin та реакція на нього');

xlabel('Час');

ylabel('Амплытуда');

 

y2=filter(b2,a2,fun);

figure(2);

plot(t,y2);hold on;

plot(t,fun);

title('Графік sin та реакція на нього');

xlabel('Час');

ylabel('Амплытуда');

 figure(3)

plot(t1,y1(1:length(t1))),hold on;

plot(t1,fun(1:length(t1)),'g');

title('Графік sin та реакція на нього');

xlabel('Час');

ylabel('Амплытуда');

 figure(4)

plot(t1,y2(1:length(t1))),hold on;

plot(t1,fun(1:length(t1)),'g');

title('Графік sin та реакція на нього');

xlabel('Час');

label('Амплытуда');

 

 
 

 

fs=256;

t=0:1/fs:0.1;

fun=sin(10*t);

ran1=rand(1,(max(length(a1),length(b1))-1));

y1=filter(b1,a1,fun,ran1);

plot(t,y1,'r')

hold on;

plot(t,fun);

ran2=rand(1,(max(length(a2),length(b2))-1));

y2=filter(b2,a2,fun,ran2);

figure(2);

plot(t,y2,'r');

hold on;

plot(t,fun);

 

 

3. Програмно визначити коефіцієнт  передачі напруги систем 1 і 2 на  частоті 10 Гц, а також різницю  фаз між вихідним і вхідним  сигналом.

fs=256;

t=0:1/fs:2;

fun=sin(10*t);

y1=filter(b1,a1,fun);

y2=filter(b2,a2,fun);

A1=max(fun);

A2=max(y1);

A3=max(y2);

h1=A2/A1;

h2=A3/A1;

 x1=fun/A1;

s1=asin(x1);

 x2=y1/A2;

s2=asin(x2);

fin1=s2-s1;

 x3=y2/A3;

j=find(asin(x3)<0);

s3=asin(x3);

fin2=s3-s1;

disp(['Коефіцієнт передачі напруги 1=',num2str(h1)]);

disp(['Коефіцієнт передачі напруги 2=',num2str(h2)]);

disp(['Різниця фаз 1=',num2str(fin1)]);

disp(['Різниця фаз 2=',num2str(fin2)]);

 

Коефіцієнт передачі напруги 1=0.13203

Коефіцієнт передачі напруги 2=0.082503

Різниця фаз 1=0+0i               0.23039+0i               0.36468+0i               0.35743+0i

Різниця фаз 2=0     0.17559     0.30606     0.33013     0.31292     0.37255

4. Розрахувати за допомогою функції filterперші 50 відліків імпульсної характеристики систем (1) та (2) подавши на вхід системи потрібний тестовий сигнал. Побудувати графіки вхідних та вихідних сигналів (функція stem), зробити висновки.

fun= [1 zeros(1, 49)];

y1=filter(b1,a1,fun);

y2=filter(b2,a2,fun);

stem([0:1:49],fun);hold on;

stem([0:1:49],y1,'r');

title('Імпульсна харак. і вхідний сигнал');

xlabel('Відліки');

ylabel('Амплітуда');

figure(2)

stem([0:1:49],fun); hold on;

stem([0:1:49],y2,'r');

title('Імпульсна харак. і вхідний сигнал');

xlabel('Відліки');

ylabel('Амплітуда');

 

 

5. Розрахувати імпульсну характеристику  систем (1) та (2) по коефіцієнтам різницевих  рівнянь з використанням функції  impz. Розрахувати 50 та 500 відліків. Порівняти результати з результатами п. 3, побудувати графіки, зробити висновки.

h1_i50=impz(b1, a1, 50);

h1_i500=impz(b1, a1, 500);

h2_i50=impz(b2, a2, 50);

h2_i500=impz(b2, a2, 500);

figure

subplot(4, 1, 1)

stem(h1_i50, '.')

title('Імпульсна характеристика, 50 відліків')

subplot(4, 1, 2)

stem(h1_i500, '.')

title('Імпульсна характеристика, 500 відліків')

subplot(4, 1, 3)

stem(h2_i50, '.')

title('Імпульсна характеристика, 50 відліків')

subplot(4, 1, 4)

stem(h2_i500, '.')

title('Імпульсна характеристика, 500 відліків')

 

6. Сформувати два синусоїдальних сигнали частоти 3 та 20 Гц тривалістю 1 с. Проілюструвати властивість адитивності системи, визначивши реакцію системи спочатку на кожний з сигналів окремо, а потім на суму цих сигналів.

Проілюструвати властивість однорідності системи.

fs=256;

t=0:1/fs:1;

fun1=sin(20*t);

fun2=sin(3*t);

y1=filter(b1,a1,fun1);

y2=filter(b1,a1,fun2);

fun3=fun1+fun2;

y3=filter(b1,a1,fun3);

plot(t,y1,'r');

hold on;

plot(t,y2);

hold on;

plot(t,y3,'g');

hold on;

title('Реакція системи(1) на синус сигнал');

xlabel('Час');

ylabel('Амплітуда');

   s=y1+y2;

  y4=y3-s;

figure(2)

plot(t,y4,'g')

title('Адетивність_1');

xlabel('Час');

ylabel('Амплітуда');

ys1=filter(b2,a2,fun1);

ys2=filter(b2,a2,fun2);

  ys3=filter(b2,a2,fun3);

figure(3)

plot(t,ys1,'r');hold on;

plot(t,ys2);hold on;

plot(t,ys3,'g');hold on;

title('Реакція системи(2) на синус сигнал');

xlabel('Час');

ylabel('Амплітуда');

   ss=ys1+ys2;

  ys4=ys3-ss;

figure(4)

plot(t,ys4,'g')

title('Адетивність_2');

xlabel('Час');

ylabel('Амплітуда');

c=5;

funo=c*fun1;

f1=filter(b1,a1,funo);

f2=filter(b2,a2,funo);

  g1=c*filter(b1,a1,fun1);

  g2=c*filter(b2,a2,fun1);

   d1=f1-g1;

  d2=f2-g2;

figure(5)

plot(t,d1)

title('Однорідність_1');

xlabel('Час');

ylabel('Амплітуда');

figure(6)

plot(t,d2)

title('Однорідність_2');

xlabel('Час');

 ylabel('Амплітуда');

 

7. Розрахувати реакцію систем (1) та (2) на сигнал з п. 2 з використанням  функції розрахунку згортки conv. Побудувати графіки вхідного та вихідного сигналу, аналогічні п. 2 (з нульовими початковими умовами). Порівняти з результатами виконання п. 2. Зробити висновки.

fs=256;

t=0:1/fs:2;

fun=sin(10*t);

[h1,t1]=impz(b1,a1);

w1=conv(fun,h1);

plot(t,w1(1:length(t)));

title('Реакція системи 2 sin на сигнал');

xlabel('Час');

ylabel('Амплітуда');

[h2,t2]=impz(b2,a2);

w2=conv(fun,h2);

figure(2)

plot(t,w2(1:length(t)));

title('Реакція системи 2 sin на сигнал');

xlabel('Час');

ylabel('Амплітуда');

y1=filter(b1,a1,fun);

y2=filter(b2,a2,fun);

s1=w1(1:length(t))-y1;

s2=w2(1:length(t))-y2;

figure(3)

plot(t,s1);

title('Графік порівняння 1');

xlabel('Час');

ylabel('Амплітуда');

 figure(4)

plot(t,s2);

title('Графік порівняння 2');

xlabel('Час');

ylabel('Амплітуда');

8. Відновити імпульсну характеристику  ЛДС (1) та (2) по відомим реакції  та вхідному сигналу з використанням  результатів п. 6 (функція deconv). Зробити висновки.

9. Обчислити комплексну частотну  характеристику систем (1) та (2) з  використанням функції freqz, побудувати графіки АЧХ та ФЧХ. Розрахувати 100 значенньКЧХ:

9.1. для частоти дискретизації  256 Гц;

9.2. для частоти дискретизації  1024 Гц.

Зробити висновки.

f = 0:500/256:500*(1-1/256);

[H1_20, w] = freqz(b1, a1, 1024);

H1_256 = freqz(b1, a1, f, 256);

[H2_20, w] = freqz(b2, a2, 1024);

H2_256 = freqz(b2, a2, f, 256);

 

freqz(b1,a1,100,256);

figure(2)

freqz(b2,a2,100,256);

figure(3)

freqz(b1,a1,100,1024);

figure(4)

freqz(b2,a2,100,1024);

 

  

10. Розрахувати АЧХ та ФЧХ  систем (1) та (2) по обчисленій в  п. 9 комплексній частотній характеристиці  з використанням функцій abs та phase. Побудуватиграфіки, порівняти з результатами п.9. Зробити висновки щодо характеру зміни модуля коефіцієнта передачи системи з частотою. Порівняти дві системи.

f = 0:500/256:500*(1-1/256);

[H1_20, w] = freqz(b1, a1, 1024);

H1_256 = freqz(b1, a1, f, 256);

[H2_20, w] = freqz(b2, a2, 1024);

H2_256 = freqz(b2, a2, f, 256);

AH1_20=abs(H1_20);

AH2_20=abs(H2_20);

AH1_256=abs(H1_256);

AH2_256=abs(H2_256);

PH1_20=phase(H1_20);

PH2_20=phase(H2_20);

PH1_256=phase(H1_256);

PH2_256=phase(H2_256);

figure

subplot(2,2,1)

plot(w, AH1_20)

subplot(2,2,2)

plot(w, AH2_20)

subplot(2,2,3)

plot(w, PH1_20)

subplot(2,2,4)

plot(w, PH2_20)

figure

subplot(2,2,1)

plot(f, AH1_256)

subplot(2,2,2)

plot(f, AH2_256)

subplot(2,2,3)

plot(f, PH1_256)

subplot(2,2,4)

plot(f, PH2_256)

 

11. Побудувати функцію, яка визначає  значення АЧХ та ФЧХ систем 1 та 2 на довільній частоті. Перевірити  за допомогою отриманої функції  правильність розрахунків з п. 3, зробити висновки.

12. Розрахувати реакцію ЛДС (1) та (2) на послідовність прямокутних  імпульсів зі шпаруватістю 30 %. Побудувати  графіки вхідного та вихідного  сигналів, зробити висновки щодо  спотворення вихідного сигналу  відносно вхідного.

t=0:1/256:2;

d=0:1/4:2;

fun=pulstran(t,d,'rectpuls',0.3);

plot(t,fun);hold on;

title('Прямокутні імпульси');

xlabel('Час');

ylabel('Амплітуда');

 

y1=filter(b1,a1,fun);

y2=filter(b2,a2,fun);

plot(t,y1,'g');

title('Реакція 1 на послідовність прям.імпульси');

xlabel('Час');

ylabel('Амплітуда');

 

figure(2)

plot(t,y2,'g');hold on;

plot(t,fun);

title('Реакція 2 на послідовність прям.імпульси');

xlabel('Час'); ylabel('Амплітуда');

 

 

13. Розрахувати реакцію ЛДС (1) та (2) на вхідний сигнал виміряної раніше власної ЕКГ тривалістю 5 секунд. Побудувати графіки вхідного та вихідного сигналів.

14. Розрахувати реакцію ЛДС на записаний звуковий сигнал. Прослухати вихідний сигнал. Зробити висновки щодо спотворень, які вносять до сигналів системи (1) та (2).

[fun,fs,nbits]=wavread('oksana.wav');

  t=0:1/fs:1;

y1=filter(b1,a1,fun);

y2=filter(b2,a2,fun);

plot(t,fun(1:length(t)));hold on;

plot(t,y1(1:length(t)),'g');

title('Звуковий сигнал 1');

xlabel('Час');

ylabel('Амплітуда');

 

figure(2)

plot(t,y2(1:length(t)),'g');hold on;

plot(t,fun(1:length(t)));

title('Звуковий сигнал 2');

xlabel('Час');

ylabel('Амплітуда');

wavplay(y1,fs);

wavplay(y2,fs);