bwthui.m
上传用户:loeagle
上传日期:2013-03-02
资源大小:1236k
文件大小:7k
源码类别:

通讯编程文档

开发平台:

Matlab

  1. function bwthui(cfreq,action)
  2. %impulse response and frequency response of Butterworth Lowpass
  3. if nargin == 0 cfreq = 5/(4*pi); end;
  4. if isstr(cfreq) cfreq = str2num(alpha); end;
  5. if (cfreq>1) cfreq = 1; end;
  6. if (cfreq<0) cfreq = 0; end;
  7. if (nargin<2)
  8.   action = 'start';
  9. end;
  10. if (strcmp(action,'start') | strcmp(action,'update'))
  11. % define constants 
  12. laenge = 5;
  13. aufloesung = 20; % Samples per unit
  14. % identical values as in  Simulink
  15. order = 4;
  16. cf = 4*tan(pi*cfreq/(aufloesung/4)/2);
  17.  
  18.     
  19.     % get zeros, poles and gain
  20.         
  21.     z = [];
  22.     p = exp(sqrt(-1)*(pi*(1:2:2*order-1)/(2*order)+pi/2)).';
  23.     k = real(prod(-p));
  24.     
  25.     % Compute transfer function
  26.     p1 = cfreq * p; 
  27.     f1=linspace(-1,1);
  28.     w1 = i*f1;
  29.     w1=repmat(w1,order,1);
  30.     p1=repmat(p1,1,length(f1));
  31.     prod(w1-p1);
  32.     Sf1 = k./prod(w1-p1);
  33.     h1 = abs(Sf1);
  34.     h1 = h1/max(h1);
  35.     
  36.     % Compute impulse response
  37.     % first transform from s-space to state space
  38.     p = p(isfinite(p));
  39.     z = z(isfinite(z));
  40.     np = length(p);
  41.     nz = length(z);
  42.     z = cplxpair(z,1e6*nz*norm(z)*eps + eps);
  43.     p = cplxpair(p,1e6*np*norm(p)*eps + eps);
  44.     a=[]; b=zeros(0,1); c=ones(1,0); d=1;
  45.     if rem(np,2) & rem(nz,2)
  46.         a = p(np);
  47.         b = 1;
  48.         c = p(np) - z(nz);
  49.         d = 1;
  50.         np = np - 1;
  51.         nz = nz - 1;
  52.     end
  53.  
  54.     if rem(np,2)
  55.         a = p(np);
  56.         b = 1;
  57.         c = 1;
  58.         d = 0;
  59.         np = np - 1;
  60.     end 
  61.     if rem(nz,2)
  62.         num = real(poly(z(nz)));
  63.         den = real(poly(p(np-1:np)));
  64.         wn = sqrt(prod(abs(p(np-1:np))));
  65.         if wn == 0, wn = 1; end
  66.             t_temp = diag([1 1/wn]); 
  67.             a = t_temp[-den(2) -den(3); 1 0]*t_temp;
  68.             b = t_temp[1; 0];
  69.             c = [1 num(2)]*t_temp;
  70.             d = 0;
  71.             nz = nz - 1;
  72.             np = np - 2;
  73.         end
  74.         v = 1;
  75.         
  76.         while v < nz
  77.             index = v:v+1;
  78.             num = real(poly(z(index)));
  79.             den = real(poly(p(index)));
  80.             wn = sqrt(prod(abs(p(index))));
  81.     
  82.             if wn == 0, wn = 1; end
  83.                 
  84.             t_temp = diag([1 1/wn]); 
  85.             a1 = t_temp[-den(2) -den(3); 1 0]*t_temp;
  86.             b1 = t_temp[1; 0];
  87.             c1 = [num(2)-den(2) num(3)-den(3)]*t_temp;
  88.             d1 = 1;
  89.             [ma1,na1] = size(a);
  90.             [ma2,na2] = size(a1);
  91.             a = [a zeros(ma1,na2); b1*c a1];
  92.             b = [b; b1*d];
  93.             c = [d1*c c1];
  94.             d = d1*d;
  95.             v = v + 2;
  96.         end
  97.         while v < np
  98.             den = real(poly(p(v:v+1)));
  99.             wn = sqrt(prod(abs(p(v:v+1))));
  100.     
  101.             if wn == 0, wn = 1; end
  102.             
  103.             t_temp = diag([1 1/wn]); 
  104.             a1 = t_temp[-den(2) -den(3); 1 0]*t_temp;
  105.             b1 = t_temp[1; 0];
  106.             c1 = [0 1]*t_temp;
  107.             d1 = 0;
  108.  
  109.             [ma1,na1] = size(a);
  110.             [ma2,na2] = size(a1);
  111.             a = [a zeros(ma1,na2); b1*c a1];
  112.             b = [b; b1*d];
  113.             c = [d1*c c1];
  114.             d = d1*d;
  115.             v = v + 2;
  116.         end
  117.     % Consider cutoff frequency
  118.     a = a * cf;
  119.     b= b * cf;
  120.  
  121.    t_temp = 1/2;
  122. r_temp = sqrt(t_temp);
  123. t1_temp = eye(size(a)) + a*t_temp/2;
  124. t2_temp = eye(size(a)) - a*t_temp/2;
  125. a = t2_tempt1_temp;
  126.     b = t_temp/r_temp*(t2_tempb);
  127.     a = poly(a);
  128.   
  129.     cf = 2*atan2(cf,4);
  130.     
  131.     r = -ones(order,1);
  132.     w = 0;
  133.     
  134.     b = poly(r);
  135.     
  136.     ehochj = exp(-j*w*(0:length(b)-1));
  137.     b = real(b*(ehochj*a(:))/(ehochj*b(:)));
  138.     
  139.     % compute impulse response
  140. t = (0:(laenge*aufloesung-1))';
  141.     imp = filter(b,a,t==0);
  142.     t = t/aufloesung;
  143.     
  144. % response to rectangular pulse
  145.     
  146. x = [ones(aufloesung,1); zeros((laenge-1)*aufloesung,1)];
  147. rec = filter(b,a,x);
  148. end;
  149. if strcmp(action,'start')
  150. % open window
  151. set(0,'Units','pixels');
  152. scnsize = get(0,'ScreenSize');
  153. figure ('Position', [0.05*scnsize(3)   0.3*scnsize(4)   0.9*scnsize(3)   0.4*scnsize(4)], ...
  154. 'Name', 'Impulse Response, Rectangular Response and Frequency Response of 4th Order Butterworth Lowpass', ...
  155. 'Tag', 'Butterworth', ...
  156. 'NumberTitle', 'off' ...
  157. );
  158. % ----------------------------------
  159. % Slider for cuttoff frequency
  160. % ----------------------------------
  161. text = uicontrol(gcf, ...
  162. 'Tag', 'BwthTextfeld', ...
  163. 'Style', 'text', ...
  164.         'Units', 'normalized', ...
  165. 'Position', [.01 .82 .08 .13], ...
  166. 'BackgroundColor', 'red', ...
  167. 'ForegroundColor', 'white', ...
  168.         'String', sprintf('Cut-offnfrequencyn%1.3f',cfreq) ...
  169. );
  170.     cb = 'bwthui(get(findobj(gcf,''Tag'',''BwthSlider''),''Value''),''update'');';
  171. slider = uicontrol(gcf, ...
  172. 'Tag', 'BwthSlider', ...
  173. 'Style', 'slider', ...
  174.         'Units', 'normalized', ...
  175. 'Position', [.04 .2 .02 .6], ...
  176. 'Min', 0, ...
  177. 'Max', 1, ...
  178. 'Value', cfreq, ...
  179. 'Callback', cb ...
  180. );
  181. % -------------------------------------------
  182. % Plot 1: impulse response
  183. % -------------------------------------------
  184. subplot(1,3,1);
  185. plot(t,imp,'EraseMode','background');
  186. set(gca, 'Tag', 'BwthImpulsantwort');
  187. set(gcf,'DefaultTextColor','m')
  188. xlabel('t/T');
  189. ylabel('g(t/T)');
  190. title('Impulse Response');
  191. set(gca,'XLimMode','manual');
  192. set(gca,'XLim', [0 laenge]);
  193. set(gca,'XTick', [0 : 1 : laenge]);
  194. set(gca,'YLimMode','manual');
  195. set(gca,'YLim', [-.025 .125]);
  196. grid
  197. % -------------------------------------------
  198. % Plot 2: response to rectangular pulse
  199. % -------------------------------------------
  200. subplot(1,3,2);
  201. plot(t,rec,t,x,'EraseMode','background');
  202. set(gca, 'Tag', 'BwthRechteckantwort');
  203. set(gcf,'DefaultTextColor','m')
  204. xlabel('t/T');
  205. ylabel('g_r(t/T)');
  206. title('Response to Rectangular Pulse');
  207. set(gca,'XLimMode','manual');
  208. set(gca,'XLim', [0 laenge]);
  209. set(gca,'XTick', [0 : 1 : laenge]);
  210. set(gca,'YLimMode','manual');
  211. set(gca,'YLim', [-.25 1.25]);
  212. grid
  213. % -----------------------
  214. % Plot 3: Frequency Response
  215. % -----------------------
  216. subplot(1,3,3);
  217. plot(f1,h1,'EraseMode','background');
  218.     
  219. set(gca, 'Tag', 'BwthFrequenzgang');
  220. set(gcf,'DefaultTextColor','m')
  221. xlabel('f T');
  222. ylabel('G(f/T)/T');
  223. title('Fourier Transform');
  224. set(gca,'XLimMode','manual');
  225. set(gca,'XLim', [-1 1]);
  226. set(gca,'YLimMode','manual');
  227. set(gca,'YLim', [ 0 1.125]);
  228. grid
  229. drawnow;
  230. end;
  231. if strcmp(action, 'update')
  232.         set(0,'CurrentFigure',findobj(0,'Tag', 'Butterworth'));
  233.         set(findobj(gcf,'Tag','BwthTextfeld'),'String', sprintf('Cut-offnfrequencyn%1.3f',cfreq));
  234. set(get(findobj(gcf,'Tag','BwthImpulsantwort'),'Children'),'YData',imp);
  235. chld = get(findobj(gcf,'Tag','BwthRechteckantwort'),'Children');
  236. set(chld(1),'YData',x); set(chld(2),'YData',rec); clear chld;
  237.     set(get(findobj(gcf,'Tag','BwthFrequenzgang'),'Children'),'YData',h1);
  238.     
  239. drawnow;
  240. end;
  241. clear action;