Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- clear all
- close all
- clc
- Fs = 6100;
- Fn = Fs / 2;
- Rp = 3;
- Rs = 25;
- Wp = [500 2000]/Fn;
- Ws = [707 1414]/Fn;
- K = 1024;
- f1=Fs/pi*tan((pi*500)/Fs)
- f2=Fs/pi*tan((pi*2000)/Fs)
- f3=Fs/pi*tan((pi*707)/Fs)
- f4=Fs/pi*tan((pi*1414)/Fs)
- figure(1)
- hold on;
- plot([0 500],[25 25]);
- plot([500 500],[0 25]);
- plot([707 1414],[3 3]);
- plot([2000 2000],[0 25]);
- plot([2000 Fn],[25 25]);
- axis([0 Fn 0 45]);
- xlabel('f,Hz');
- ylabel('a,dB');
- title('Габарит на ЛФ');
- [N,Wn] = buttord(Wp,Ws,Rp,Rs);
- [Nz,Dz] = butter(N,Wn);
- figure(2)
- hold on;
- plot([0 500],[25 25]);
- plot([500 500],[0 25]);
- plot([707 1414],[3 3]);
- plot([2000 2000],[0 25]);
- plot([2000 Fn],[25 25]);
- axis([0 Fn 0 45]);
- [H,w] = freqz(Nz,Dz);
- ma = -20*log10(abs(H));
- plot(w*Fs/(2*pi),ma,'r-.');
- xlabel('f,Hz');
- ylabel('a,dB');
- title('Габарит на ЛВ и крива на затихване');
- legend('Габарит','Затихване');
- figure(3)
- [z,p,k] = butter(N,Wn);
- subplot(211);
- zplane(Nz,Dz);
- title('ПНД');
- fi = angle(H);
- fi1 = unwrap(fi);
- subplot(212)
- plot(w*Fs/(2*pi), fi1*180/pi);
- title('ФЧХ');
- xlabel('Честота (rad/s)');
- ylabel('ФЧХ (rad)');
- figure(4)
- [h,t] = impz(Nz,Dz);
- stem(t,h,'m');
- grid;
- title('Импулсна характеристика');
- xlabel('Отчети на времето');
- ylabel('Амплитуда (volts)');
- figure(5)
- n = (0:199)/Fs;
- x = sin(2*pi*1000*n);
- y = filter(Nz,Dz,x);
- subplot(221);
- plot(n,x);
- title('Входен сигнал');
- xlabel('Амплитуда (volts)');
- ylabel('Отчети на времето');
- subplot(223);
- plot(n,y);
- title('Изходен сигнал');
- xlabel('Амплитуда (volts)');
- ylabel('Отчети на времето');
- Px = fft(x,K);
- px = abs(Px(1:(K/2)));
- Py = fft(y,K);
- py = abs(Py(1:(K/2)));
- subplot(222);
- f = (0:(length(Px) - 1) / 2);
- plot(f,px);
- xlabel('Честота (Hz)');
- title('Спектър на входен сигнал');
- subplot(224)
- f1 = (0:(length(Py)-1)/2);
- plot(f1,py);
- xlabel('Честота (Hz)');
- title('Спектър на изходен сигнал');
- b = [0.0177 0 -0.0887 0 0.1774 0 -0.1774 0 0.0887 0 -0.0177];
- a = [1.0000 -3.0507 5.0514 -5.9813 5.7213 -4.3294 2.5774 -1.1966 0.4323 -0.1060 0.0151];
- [b0,B,A] = dir2cas(b,a)
- figure(6)
- n1 = (0:199)/Fs;
- x1 = sin(2*pi*250*n1);
- y1 = filter(Nz,Dz,x1);
- subplot(221);
- plot(n1,x1);
- title('Входен сигнал');
- xlabel('Амплитуда (volts)');
- ylabel('Отчети на времето');
- subplot(223);
- plot(n1,y1);
- title('Изходен сигнал');
- xlabel('Амплитуда (volts)');
- ylabel('Отчети на времето');
- Px1 = fft(x1,K);
- px1 = abs(Px1(1:(K/2)));
- Py1 = fft(y1,K);
- py1 = abs(Py1(1:(K/2)));
- subplot(222);
- f2 = (0:(length(Px1) - 1) / 2);
- plot(f2,px1);
- xlabel('Честота (Hz)');
- title('Спектър на входен сигнал');
- subplot(224)
- f3 = (0:(length(Py1)-1)/2);
- plot(f3,py1);
- xlabel('Честота (Hz)');
- title('Спектър на изходен сигнал');
- figure(7)
- hold on;
- %WsNew=373.92;
- %WpNew=3450.15;
- Wp1 = [511 3233];
- Ws1 = [740 1731];
- [N2,Wn2]=buttord(Wp1,Ws1,Rp,Rs,'s')
- [Ns,Ds]=butter(N2,Wn2,'s')
- printsys(Ns,Ds)%принтира Ns и Ds
- [Ts,ws]=freqs(Ns,Ds);%изчислява аналоговата предавателна функция
- ma1=-20*log10(abs(Ts));%затихването в децибели
- plot(ws,ma1,'g--'); %крива на затихване на аналогов ЛФ
- plot(w*Fs/(2*pi),ma,'r-.');%крива на затихване на цифров ЛФ
- ylabel('a,dB');
- xlabel('f,Hz');
- title('Криви на затихването на аналогов ЛФ и цифров ЛФ');
- legend('Крива на затихване на аналогов ЛФ','Крива на затихване на цифров ЛФ');
- plot([0 500],[25 25]);
- plot([500 500],[0 25]);
- plot([707 1414],[3 3]);
- plot([2000 2000],[0 25]);
- plot([2000 Fn],[25 25]);
- axis([0 Fn 0 45]);
- xlabel('f,Hz');
- ylabel('a,dB');
- title('Габарит на ЛФ');
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement