Статьи о MathLab / LAB1EXMP
.DOCЦифровые фильтры и их свойства
Курсовая работа состоит из теоретической и практической частей, выводов по работе и приложения. Теоретическая часть работы содержит изложение основ применения дискретного преобразования Фурье для спектрального анализа сигналов и общие основы теории цифровой фильтрации с подробным описанием метода проектирования цифрового КИХ-фильтра методом окна. Практическая часть работы представлена функциями 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, построить график модуля спектральной функции.
Использовать окна, задаваемые функцией:
-
BARTLETT(N) returns the N-point Bartlett window.
-
BLACKMAN(N) returns the N-point Bartlett window.
-
BOXCAR(N) returns the N-point rectangular window.
-
CHEBWIN(N,R) returns the N-point Chebyshev window with R decibels of ripple.
-
HAMMING(N) returns the N-point Hamming window.
-
HANNING(N) returns the N-point Hanning window.
-
KAISER(N,beta) returns the BETA-valued N-point Kaiser window.
-
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))]);