ID: 51128

Название работы: Фільтрація сигналів

Категория: Лабораторная работа

Предметная область: Информатика, кибернетика и программирование

Описание: Мета роботи: набути навичок проектування цифрових фільтрів, задавання специфікації фільтрів залежно від властивостей сигналів, які треба фільтрувати; набути навичок реалізації дискретної фільтрації сигналів у середовищі MatLAB.

Язык: Украинкский

Дата добавления: 2014-02-06

Размер файла: 889.81 KB

Работу скачали: 17 чел.

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

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

Факультет електроніки

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

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

«Фільтрація сигналів»

Виконав:  студент 3-го курсу

групи ДП-92

 Лонтковський С.А.

Київ – 2011

Мета роботи: набути навичок проектування цифрових фільтрів, задавання специфікації фільтрів залежно від властивостей сигналів, які треба фільтрувати; набути навичок реалізації дискретної фільтрації сигналів у середовищі MatLAB.

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

1. Сформувати вектор відліків часу тривалістю 5 с для частоти дискретизації 128 Гц. Сформувати прямокутний імпульс в момент часу 3 с тривалості 0.1 с. Спроектувати ФНЧ Батерворта для позбавлення сигналу від шуму (функції buttord, butter, filter). Оцінити відношення сигнал/шум на вході та виході фільтра.

2. Сформувати вектор відліків часу тривалістю 1 с для частоти дискретизації 128 Гц. Сформувати сигнали ділянки синусоїди частотою 10 Гц амплітуди 1 В. Додати випадковий сигнал з нульовим середнім значенням амплітуди 2 В. Спроектувати ФНЧ Чебишова І-го роду, ФВЧ Чебишова ІІ-роду, СФ Кауера для позбавлення сигналу від шуму. Оцінити відношення сигнал/шум на вході та виході фільтра.

4. Для оцифрованих сигналів електрокардіограми, електроенцефалограми, прочитаної з файлу, а також ЕЕГ здорової та хворої людини, сигналів артеріального та внутрішньочерепного тиску та плетизмограми спроектувати фільтри для позбавлення від мережевої перешкоди 50 Гц.

5. Для звукових сигналів, які отримані з різною частотою дискретизації, виконати розділення на три спектральні діапазони: до 450 Гц; від 450 Гц до 1 кГц; від 1 кГц до 4 кГц з використанням фільтрів Кауера. Прослухати отримані сигнали, зробити висновки.

close all;

clear all;

clc;

 

%завдання 1------------------------------------------------

fs=128;

T=5;

fn=fs/2;

Wp=15/fn;Ws=20/fn;

Rp=3;

Rs=25;

f=[0:1/T:fs-1/T];

t=[0:1/fs:T-1/fs];

x=rectpuls(t-3,0.1);

ssx=abs(fft(x)/length(x));

figure()

subplot(2,2,1);plot(t,x);grid on;

title('Вхідний сигнал X');xlabel('Час, c');ylabel('Значення');

subplot(2,2,2);plot(f,ssx);grid on;

title('Амплітудний спектр сигнала x');

xlabel('Частота, Гц');ylabel('Значення');

N=rand(size(x))-0.5;

subplot(2,2,3);plot(t,N);grid on;title('Шум ');

xlabel('Час, c');ylabel('Значення');

ssN=abs(fft(N))/length(N);

subplot(2,2,4);plot(f,ssN);grid on;

title('Амплітудний спектр шума');

xlabel('Частота, Гц');ylabel('Значення');

Nx=x+N;

figure()

subplot(2,2,1);plot(t,Nx);grid on;title('Зашумленний вхідний сигнал');

xlabel('Час, с');ylabel('Значення');

ssNx=abs(fft(Nx)/length(Nx));

subplot(2,2,2);plot(f,ssNx);grid on;title('Амплітудний спектр Nx');

xlabel('Частота, Гц');ylabel('Значення');

[N,W]=buttord(Wp,Ws,Rp,Rs);

[b,a]=butter(N,W);

h=freqz(b,a,f,fs);

subplot(2,1,2);plot(f,abs(h));grid on;title('Амплітудний спектр Фільтра');

xlabel('Частота, Гц');ylabel('Значення');

y=filter(b,a,Nx);

figure()

subplot(2,1,1);plot(t,y);grid on;title('Відфільтрований сигнал');

xlabel('Час, с');ylabel('Значення');

ssy=abs(fft(y)/length(y));

subplot(2,1,2);plot(f,ssy);xlabel('Частота, Гц');ylabel('Значення');

title('Амплітудний спектр відфільтрованого сигналу');grid on;

 

%завдання 2------------------------------------------------

fs=128;

T=1;

fn=fs/2;

Wp=15/fn;

Ws=20/fn;

Rp=3;

Rs=70;

Rs3=150;

t=[0:1/fs:T-1/fs];

f=[0:1/T:fs-1/T];

x=sin(2*pi*10*t);

figure()

subplot(2,2,1);plot(t,x);grid on;title('Вхідний синусоїдальний сигнал');

xlabel('Час, с');ylabel('Значення');

ssx=abs(fft(x)/length(x));

subplot(2,2,2);plot(t,ssx);grid on;

title('Амплітудний спектр синусоїдального сигнала');

xlabel('Частота, Гц');ylabel('Значення');

r=2*rand(1,length(t));

subplot(2,2,3);plot(t,r);grid on;title('Випадковий сигнал');

xlabel('Час, с');ylabel('Значення');

ssr=abs(fft(r)/length(r));

subplot(2,2,4);plot(f,ssr);grid on;

title('Амплітудний спектр випадкового сигналу');

xlabel('Частота, Гц');ylabel('Значення');

rx=x+r;

figure()

subplot(2,2,1);plot(t,rx);grid on;

title('Вхідний і випадковий сигнал');

xlabel('Час, с');ylabel('Значення');

ssrx=abs(fft(rx)/length(rx));

subplot(2,2,2);plot(ssrx);title('Амплітудний спектр сигнала rs');

xlabel('Частота, Гц');ylabel('Значення');grid on;

[N1,Wn1]=cheb1ord(Wp,Ws,Rp,Rs);

[b,a]=cheby1(N1,Rp,Wn1,'low');

h=freqz(b,a,f,fs);

subplot(2,1,2);plot(f,abs(h));grid on;title('Чебешов 1-го порядку');

xlabel('Частота, Гц');ylabel('Значення');

y1=filter(b,a,rx);

figure()

subplot(2,2,1);plot(t,y1);grid on;

title('Перший відфільтрований сигнал');

xlabel('Час, с');ylabel('Значення');

ssy1=abs(fft(y1)/length(y1));

subplot(2,2,2);plot(f,ssy1);grid on;

title('Амплітудний спектр');xlabel('Частота, Гц');ylabel('Значення');

[N2,Wn2]=cheb2ord(Wp,Ws,Rp,Rs);

[b2,a2]=cheby2(N2,Rp,Wn2,'high');

y2=filter(b2,a2,rx);

subplot(2,2,3);plot(t,y2);grid on;

title('Другий відфільтрований сигнал');

xlabel('Час, с');ylabel('Значення');

ssy2=abs(fft(y2)/length(y2));

subplot(2,2,4);plot(f,ssy2);grid on;

title('Амплітудний спектр');xlabel('Частота, Гц');ylabel('Значення');

h2=freqz(b2,a2,f,fs);

figure()

subplot(2,2,1);plot(f,abs(h2));grid on;title('Чебешов 2-го порядка');

xlabel('Частота, Гц');ylabel('Значення');

[N3,Wn3]=ellipord([15/fn 17/fn],[13/fn 20/fn],Rp,Rs3);

[b3,a3]=ellip(N3,Rp,Rs,Wn3);

h3=freqz(b3,a3,f,fs);

subplot(2,2,2);plot(f,abs(h3));title('СФ Кауера');grid on;

xlabel('Частота, Гц');ylabel('Значення');

y3=filter(b3,a3,rx);

subplot(2,2,3);plot(t,y3);title('Третій відфільтрований сигнал');

xlabel('Час, с');ylabel('Значення');grid on;

ssy3=abs(fft(y3)/length(y3));

subplot(2,2,4);plot(ssy3);grid on;

title('Амплiтудний спектр');

xlabel('Частота, Гц');ylabel('Значення');

 

%завдання 4------------------------------------------------

load('ecg_10')

d=d/max(abs(d));

fs=400;

fn=fs/2;

t=[0:1/fs:(length(d)-1)/fs];

f=[0:1/t(length(d)):fn];

ssd=fft(d)/length(d);

assd=abs(ssd);

Ws=[49.9/fn 50.1/fn];Wp=[49/fn 51/fn];

[N,Wn]=ellipord(Wp,Ws,3,160);

[B,A]=ellip(N,0.5,30,Wn,'stop');

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,d);

ssy=fft(y)/length(y);

ay=abs(ssy);

subplot(4,1,1);plot(f,modh(1:length(f)));grid on;

xlabel('Частота,Гц');ylabel('Значення, мВ'); title('ЗФ фільтр');

subplot(3,2,3);plot(t,d);xlabel('Час,сек');ylabel('Значення, мВ'); title('ЕКГ');grid on;

subplot(3,2,4);plot(f,assd(1:length(f)));xlabel('Частота,Гц');ylabel('Значення, мВ');

title('Амплітудний спектр ЕКГ');grid on;

subplot(3,2,5);plot(t,y),xlabel('Час,сек');ylabel('Значення, мВ'); grid on;

subplot(3,2,6);plot(f,ay(1:length(f)));xlabel('Час,Гц');

ylabel('Значення, мВ'); grid on;

 

load('eeg_healthy_10')

sig=sig/max(abs(sig));

fs=256;

fn=fs/2;

t=[0:1/fs:(length(sig)-1)/fs];

f=[0:1/t(length(sig)):fn];

ssd=fft(sig)/length(sig);

assd=abs(ssd);

Ws=[49.9/fn 50.1/fn];Wp=[49/fn 51/fn];

[N,Wn]=ellipord(Wp,Ws,3,160);

[B,A]=ellip(N,0.5,30,Wn,'stop');

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,sig);

ssy=fft(y)/length(y);

ay=abs(ssy);

figure

subplot(4,1,1);plot(f,modh(1:length(f)));xlabel('Частота,Гц');

ylabel('Значення, мВ');title('ЗФ фільтр');grid on;

subplot(3,2,3);plot(t,sig);xlabel('Час,сек');ylabel('Значення, мВ');

title('ЕЕГ здорової');grid on;

subplot(3,2,4);plot(f,assd(1:length(f)));xlabel('Час,Гц');

ylabel('Значення, мВ');title('амплітудний спектр ЕЕГ');grid on;

subplot(3,2,5);plot(t,y);xlabel('Час,сек');ylabel('Значення, мВ');grid on;

subplot(3,2,6);plot(f,ay(1:length(f)));xlabel('Час,Гц');

ylabel('Значення, мВ');grid on;

 

load('eeg_sick_10')

sig=sig/max(abs(sig));

fs=256;

fn=fs/2;

t=[0:1/fs:(length(sig)-1)/fs];

f=[0:1/t(length(sig)):fn];

ssd=fft(sig)/length(sig);

assd=abs(ssd);

Ws=[49.9/fn 50.1/fn];Wp=[49/fn 51/fn];

[N,Wn]=ellipord(Wp,Ws,3,160);

[B,A]=ellip(N,0.5,30,Wn,'stop');

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,sig);

ssy=fft(y)/length(y);

ay=abs(ssy);

figure

subplot(4,1,1);plot(f,modh(1:length(f)));xlabel('Частота,Гц');

ylabel('Значення, мВ');title('ЗФ фільтр');grid on;

subplot(3,2,3);plot(t,sig);xlabel('Час,сек');ylabel('Значення, мВ');

title('ЕЕГ Хворої');grid on;

subplot(3,2,4);plot(f,assd(1:length(f)));xlabel('Час,Гц');

ylabel('Значення, мВ');title('амплітудний спектр ЕЕГ');grid on;

subplot(3,2,5);plot(t,y);xlabel('Час,сек');ylabel('Значення, мВ');grid on;

subplot(3,2,6);plot(f,ay(1:length(f)));xlabel('Час,Гц');

ylabel('Значення, мВ');grid on;

 

 

fid=fopen('TBI_ICP.txt', 'r');

d=fscanf(fid,'%f');

fs=125;

fn=fs/2;

t=[0:1/fs:(length(d)-1)/fs];

f=[0:1/t(length(d)):fn];

ssd=fft(d)/length(d);

assd=abs(ssd);

Ws=[49.9/fn 50.1/fn];Wp=[49/fn 51/fn];

[N,Wn]=ellipord(Wp,Ws,3,160);

[B,A]=ellip(N,0.5,30,Wn,'stop');

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,d);

ssy=fft(y)/length(y);

ay=abs(ssy);

subplot(4,1,1);plot(f,modh(1:length(f)));grid on;

xlabel('Частота,Гц');ylabel('Значення'); title('ЗФ фільтр');

subplot(3,2,3);plot(t,d);xlabel('Час,сек');ylabel('Значення');

title('Внутрішньочерепний тиск');grid on;

subplot(3,2,4);plot(f,assd(1:length(f)));xlabel('Частота,Гц');ylabel('Значення');

title('Амплітудний спектр');grid on;

subplot(3,2,5);plot(t,y),xlabel('Час,сек');ylabel('Значення'); grid on;

subplot(3,2,6);plot(f,ay(1:length(f)));xlabel('Час,Гц');

ylabel('Значення'); grid on;

 

 

fid=fopen('TBI_ABP.txt', 'r');

d=fscanf(fid,'%f');

fs=125;

fn=fs/2;

t=[0:1/fs:(length(d)-1)/fs];

f=[0:1/t(length(d)):fn];

ssd=fft(d)/length(d);

assd=abs(ssd);

Ws=[49.9/fn 50.1/fn];Wp=[49/fn 51/fn];

[N,Wn]=ellipord(Wp,Ws,3,160);

[B,A]=ellip(N,0.5,30,Wn,'stop');

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,d);

ssy=fft(y)/length(y);

ay=abs(ssy);

subplot(4,1,1);plot(f,modh(1:length(f)));grid on;

xlabel('Частота,Гц');ylabel('Значення'); title('ЗФ фільтр');

subplot(3,2,3);plot(t,d);xlabel('Час,сек');ylabel('Значення');

title('Артеріальний тиск');grid on;

subplot(3,2,4);plot(f,assd(1:length(f)));xlabel('Частота,Гц');ylabel('Значення');

title('Амплітудний спектр');grid on;

subplot(3,2,5);plot(t,y),xlabel('Час,сек');ylabel('Значення'); grid on;

subplot(3,2,6);plot(f,ay(1:length(f)));xlabel('Час,Гц');

ylabel('Значення'); grid on;

 

 

fid=fopen('vavreschuk','r');

k=fread(fid,[9,inf],'int16');

d=k(8,:);

fs=155.1;

fn=fs/2;

t=[0:1/fs:(length(d)-1)/fs];

f=[0:1/t(length(d)):fn];

ssd=fft(d)/length(d);

assd=abs(ssd);

Ws=[49.9/fn 50.1/fn];Wp=[49/fn 51/fn];

[N,Wn]=ellipord(Wp,Ws,3,160);

[B,A]=ellip(N,0.5,30,Wn,'stop');

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,d);

ssy=fft(y)/length(y);

ay=abs(ssy);

subplot(4,1,1);plot(f,modh(1:length(f)));grid on;

xlabel('Частота,Гц');ylabel('Значення, мВ'); title('ЗФ фільтр');

subplot(3,2,3);plot(t,d);xlabel('Час,сек');ylabel('Значення, мВ');

title('Плетизмограма');grid on;

subplot(3,2,4);plot(f,assd(1:length(f)));xlabel('Частота,Гц');ylabel('Значення, мВ');

title('Амплітудний спектр');grid on;

subplot(3,2,5);plot(t,y),xlabel('Час,сек');ylabel('Значення, мВ'); grid on;

subplot(3,2,6);plot(f,ay(1:length(f)));xlabel('Час,Гц');

ylabel('Значення, мВ'); grid on;

 

 

%завдання 5------------------------------------------------

[x,fs,bits]=wavread('8kHz.wav');

fn=fs/2;

t=[0:1/fs:(length(x)-1)/fs];

f=[0:1/t(length(x)):fn];

ssx=fft(x)/length(x);

assx=abs(ssx);

figure()

subplot(2,1,1);plot(t,x);xlabel('Час, сек');ylabel('Значення');

title('Вхідний сигнал');grid on;

subplot(2,1,2);plot(f,assx(1:length(f)));xlabel('Частота, Гц');

ylabel('Значення');title('Амплітудний спектр вхідного сигналу');grid on;

[N,W]=cheb2ord(450/fn,460/fn,5,10);

[B,A]=cheby2(N,30,W);

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,x);

ssy=fft(y)/length(y);

assy=abs(ssy);

wavplay(y,fs)

figure()

subplot(3,1,1);plot(f,modh);xlabel('Частота, Гц');ylabel('Значення');

title('АЧХ ФНЧ(до 450Гц)');grid on;

subplot(3,1,2);plot(t,y);xlabel('Час, сек');ylabel('Значення');

title('Відфільтрований сигнал');grid on;

subplot(3,1,3);plot(f,assy(1:length(f)));xlabel('Частота, Гц');grid on;

ylabel('Значення');title('Амплітудний спектр відфільтрованого сигналу');

 

[N,W]=cheb2ord([450/fn 1000/fn],[445/fn 1005/fn],5,6);

[B,A]=cheby2(N,50,W);

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,x);

ssy=fft(y)/length(y);

assy=abs(ssy);

wavplay(y,fs)

figure()

subplot(3,1,1);plot(f,modh);xlabel('Частота, Гц');ylabel('Значення');

title('АЧХ ПФ(від 450Гц до 1000Гц)');grid on;

subplot(3,1,2);plot(t,y);xlabel('Час, сек');ylabel('Значення');

title('Відфільтрований сигнал');grid on;

subplot(3,1,3);plot(f,assy(1:length(f)));xlabel('Частота, Гц');grid on;

ylabel('Значення');title('Амплітудний спектр відфільтрованого сигналу');

 

[N,W]=cheb2ord([1000/fn 3980/fn],[990/fn 3990/fn],5,6);

[B,A]=cheby2(N,40,W);

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,x);

ssy=fft(y)/length(y);

assy=abs(ssy);

wavplay(y,fs)

figure()

subplot(3,1,1);plot(f,modh);xlabel('Частота, Гц');ylabel('Значення');

title('АЧХ ПФ(від 450Гц до 1000Гц)');grid on;

subplot(3,1,2);plot(t,y);xlabel('Час, сек');ylabel('Значення');

title('Відфільтрований сигнал');grid on;

subplot(3,1,3);plot(f,assy(1:length(f)));xlabel('Частота, Гц');grid on;

ylabel('Значення');title('Амплітудний спектр відфільтрованого сигналу');

[x,fs,bits]=wavread('44.1kHz.wav');

fn=fs/2;

t=[0:1/fs:(length(x)-1)/fs];

f=[0:1/t(length(x)):fn];

ssx=fft(x)/length(x);

assx=abs(ssx);

figure()

subplot(2,1,1);plot(t,x);xlabel('Час, сек');ylabel('Значення');

title('Вхідний сигнал');grid on;

subplot(2,1,2);plot(f,assx(1:length(f)));xlabel('Частота, Гц');

ylabel('Значення');title('Амплітудний спектр вхідного сигналу');grid on;

[N,W]=cheb2ord(450/fn,460/fn,5,10);

[B,A]=cheby2(N,30,W);

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,x);

ssy=fft(y)/length(y);

assy=abs(ssy);

wavplay(y,fs)

figure()

subplot(3,1,1);plot(f,modh);xlabel('Частота, Гц');ylabel('Значення');

title('АЧХ ФНЧ(до 450Гц)');grid on;

subplot(3,1,2);plot(t,y);xlabel('Час, сек');ylabel('Значення');

title('Відфільтрований сигнал');grid on;

subplot(3,1,3);plot(f,assy(1:length(f)));xlabel('Частота, Гц');grid on;

ylabel('Значення');title('Амплітудний спектр відфільтрованого сигналу');

 

[N,W]=cheb2ord([450/fn 1000/fn],[445/fn 1005/fn],5,10);

[B,A]=cheby2(N,50,W);

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,x);

ssy=fft(y)/length(y);

assy=abs(ssy);

wavplay(y,fs)

figure()

subplot(3,1,1);plot(f,modh);xlabel('Частота, Гц');ylabel('Значення');

title('АЧХ ПФ(від 450Гц до 1000Гц)');grid on;

subplot(3,1,2);plot(t,y);xlabel('Час, сек');ylabel('Значення');

title('Відфільтрований сигнал');grid on;

subplot(3,1,3);plot(f,assy(1:length(f)));xlabel('Частота, Гц');grid on;

ylabel('Значення');title('Амплітудний спектр відфільтрованого сигналу');

 

[N,W]=cheb2ord([1000/fn 3980/fn],[990/fn 3990/fn],5,8);

[B,A]=cheby2(N,40,W);

modh=abs(freqz(B,A,length(f),fs));

y=filter(B,A,x);

ssy=fft(y)/length(y);

assy=abs(ssy);

wavplay(y,fs)

figure()

subplot(3,1,1);plot(f,modh);xlabel('Частота, Гц');ylabel('Значення');

title('АЧХ ПФ(від 450Гц до 1000Гц)');grid on;

subplot(3,1,2);plot(t,y);xlabel('Час, сек');ylabel('Значення');

title('Відфільтрований сигнал');grid on;

subplot(3,1,3);plot(f,assy(1:length(f)));xlabel('Частота, Гц');grid on;

ylabel('Значення');title('Амплітудний спектр відфільтрованого сигналу');

Графіки.

1)

2)

4)

5)

Ч    Частота 8кГц.

Ч  Частота 44.1кГц.