Добавил:
Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:
Скачиваний:
34
Добавлен:
01.05.2014
Размер:
306.18 Кб
Скачать

Цифровые фильтры и их свойства

Курсовая работа состоит из теоретической и практической частей, выводов по работе и приложения. Теоретическая часть работы содержит изложение основ применения дискретного преобразования Фурье для спектрального анализа сигналов и общие основы теории цифровой фильтрации с подробным описанием метода проектирования цифрового КИХ-фильтра методом окна. Практическая часть работы представлена функциями MATLAB выполняющих операции представленные в задании. Выводы по работе заключаются в сравнительной оценки результатов фильтрации различных сигналов. В приложении к курсовой работе приводится описание используемых функций MATLAB.

1. Генерация входной последовательности сигнала

Ввести входной сигнал и построить его график.

Частота дискретизации: 128Гц

Рабочий диапазон частот: 0.5-35Гц

Длина реализации: 512 отсчетов (4c)

Fs = 128;

T = 4; % time

F1 = 0.5;

F2 = 35;

n = 512; % length of input signal

indx = (0:n-1)'; % index vector

Примеры задания входной последовательности

1. Полигармонический сигнал

Saicos(wit), wi = 3, 7, 10, 25Гц

Пример 1.1.

2. Полигармонический сигнал и cлучайный шум (сигнал/шум = 20дБ)

Saicos(wit) +h(t), wi = 3, 7, 10, 25Гц

format short

t = 1:n;

x = 30*cos((3*2*pi/n)*t) + 60*cos((7*2*pi/n)*t) + 90*cos((10*2*pi/n)*t) + 30*cos((25*2*pi/n)*t);

Пример 1.2.

3. Модулированный полигармонический сигнал

A(1 +Saicos(wit))sinWt

Пример 1.3.

4. Сигнал составленный из подпоследовательностей различной частоты (3, 7, 10, 25Гц)

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

4.2. равное число периодов задания подпоследовательностей сигнала

4.3. равные энергетические участки задания подпоследовательностей сигнала

Пример 1.4.1.

5. Сигнал с непрерывно возрастающей частотой

5.1. от 0.5 до 35Гц

5.2. от 0.5 до 25Гц

5.3. от 3 до 35Гц

5.4. от 3 до 25Гц

Пример 1.5.1.

6. Реальный сигнал (последовательность из 2560 16-разрядных чисел со знаком, Fs=256Гц):

6.1. signal1.dig

6.2. signal2.dig

6.3. signal3.dig

6.4. signal4.dig

6.5. signal5.dig

Пример 1.6.1.

Ввод сигнала из файла и его отображение выполняется с использованием функций: FOPEN (open file), FREAD (read binary data from file), FCLOSE (Close file), CLF (clear current figure), PLOT (linear plot).

??? format compact;Â

|

Missing variable or function.

Описание используемых функций содержится в приложении.

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

%fid = fopen('c:\temp\signal1.dig', 'r');

%x = fread(fid, n, 'bit16', 16);

%st = fclose(fid);

clf, plot(indx,x)

ylabel('Amplitude (mkV)')

xlabel('Time')

title('Input sequence')

2. Реализация спектрального анализа с использованием БПФ

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

Использовать окна, задаваемые функцией:

  1. BARTLETT(N) returns the N-point Bartlett window.

  2. BLACKMAN(N) returns the N-point Bartlett window.

  3. BOXCAR(N) returns the N-point rectangular window.

  4. CHEBWIN(N,R) returns the N-point Chebyshev window with R decibels of ripple.

  5. HAMMING(N) returns the N-point Hamming window.

  6. HANNING(N) returns the N-point Hanning window.

  7. KAISER(N,beta) returns the BETA-valued N-point Kaiser window.

  8. TRIANG(N) returns the N-point triangular window.

Пример 2.1.

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

Описание функций FFT и IFFTсодержится в приложении.

% time domain

val = get(popup,'Value');

if (val == 1),

f=amp*sin(freq*t*2*pi);

elseif (val == 2), % square wave

tt=freq*t*2*pi;

tmp=rem(tt,2*pi);

f=amp*(2*rem((tt<0)+(tmp>pi | tmp<-pi)+1,2)-1);

elseif (val == 3), % sawtooth

tt=freq*t*2*pi;

f=amp*((tt < 0) + rem(tt,2*pi)/2/pi - .5)*2;

end;

% frequency domain

F=fft(window.*f,2*M);

F=F(1:M);

w=(0:M-1)*pi/M;

FF=20*log10(abs(F));

ind=find(FF<min_dB);

FF(ind)=NaN*ones(size(ind)); % put NaN's in where

% min_dB shows up - this is to work around no clipping in xor mode

freq_line=plot(w/2/pi/T,FF,'EraseMode','xor');

axis([0 1/(2*T) min_dB 50]);

grid on;

ylabel('Magnitude (dB)');

xlabel('Frequency (Hertz)');

ax_time=axes('Position',[.12 .58 .6 .3],'XLim',[t0 t1],'YLim',[-1 1]);

time_line=plot(t,f,'EraseMode','xor');

axis([t0 t1 -1 1]);

% (set to xor mode to prevent re-rendering, that is, for speed)

grid on;

ylabel('Waveform');

xlabel('Time (Seconds)');

text(0.12, 1.55,' Click and drag waveform to change');

text(0.12, 1.3,'fundamental frequency and amplitude');

set(time_line,'ButtonDownFcn','sigdemo1(''down'')');

SIGDEMO1_DAT = [freq; amp; N; M; min_dB; 0; 0; ...

time_line; freq_line; freq_field; popup; -1; gcf; t(:); window(:)];

ADDIT_DAT = winHndl;

watchoff(oldFigNumber);

format long e

F1=fft(bartlett(n).*x');

F2=fft(x);

F1=F1(1:35/(1/T));

F2=F2(1:35/(1/T));

clf, plot(indx(1:35/(1/T)),abs(F1),'r',indx(1:35/(1/T)),abs(F2),'b')

% axis([0 128 0 10000])

axis([0 35 0 max(abs(F2))]);

grid on

ylabel('Magnitude')

xlabel('Frequency (Hertz)')

title('Power spectrum')

??? Error using ==> .*

Matrix dimensions must agree.

[b,a]=cheby1(5,0.5,pi*13/128);

[H,w]=freqz(b,a,n);

clf,plot (w,abs(H))

%axis ([0 0.5 0 max(abs(H))])

Фильтры Баттерворта

Фильтры Баттерворта нижних частот характеризуются тем, что имеют максимально гладкую амплитудную характеристику в начале координат в S-плоскости. Квадрат амплитудной характеристики нормированного (т.е. имеющего частоту среза 1 рад/с) фильтра Баттерворта равен:

,

где n - порядок фильтра.

Можно сформулировать несколько свойств фильтров Баттерворта нижних частот.

· Фильтры Баттерворта имеют только полюсы (все нули передаточных функций этих фильтров расположены на бесконечности).

· На частоте W=1 рад/с коэффициент передачи фильтра равен (т.е. на частоте среза их амплитудная характеристика спадает на 3 дБ).

· Порядок фильтра n полностью определяет весь фильтр.

[b,a] = butter(5,.5)

b contains the numerator coefficients, a contains the denominator coefficients.

The filter's Z-transform can be computed around the unit circle with the freqz function

[H,w] = freqz(b,a);

plot(w,abs(H))

ylabel('Magnitude')

xlabel('Frequency (radians)')

title('Lowpass Filter Response')

format long e

y = filter(b,a,x);

clf, plot(t,x,t,y)

wi=bartlett(n);

w=fft(wi);

w=w(1:10/(1/T));

clf,plot(abs(w))

axis([0 10 0 max(abs(w))]);