[已解决] 错误使用 symengine>@()0.0,具体是EigenFilter的Equiripple设计问题

[复制链接]
Alison_Ni 发表于 2022-6-14 20:31:30
本帖最后由 Alison_Ni 于 2022-6-14 20:40 编辑

之前是看了下面这个帖子解决的函数矩阵的问题,但是现在不能对矩阵的一个函数元素做积分运算,不知道是怎么回事。代码部分只要改odd部分就可以,even部分我后面自己再改。
https://www.ilovematlab.cn/thread-606671-1-1.html

错误信息:
  1. 错误使用 symengine>@()0.0
  2. 输入参数太多。
  3. 出错 quad (line 67)
  4. y = f(x, varargin{:});
  5. 出错 EigenFilter_Equiripple (line 66)
  6.                 P(n,m) = (quad(Epf,0,wp,0.000001) + quad(Esf,ws ,pi,0.000001 )) ;

  7. 具体出错地方是 P(n,m) = (quad(Epf,0,wp,0.000001) + quad(Esf,ws ,pi,0.000001 )) ;
复制代码
以下是原代码:
  1. %% Eigen Filter Design
  2. % example 6: lowpass EigenFilter with equiripple   

  3. clear all;clc;close all;

  4. N = 29;                     % odd
  5. fp = 0.15;
  6. fs = 0.2;
  7. alph = 0.5;
  8. wp = fp*2*pi;
  9. ws = fs*2*pi;
  10. % f = linspace(0,0.5,1001);
  11. mode = mod(N,2);            % choice odd or even mode

  12. % need revising
  13. if(mode)
  14.     M = (N-1)/2;            % length is odd
  15. else
  16.     M = N/2;                % length is even
  17. end

  18. c = ones(M+1,1);            % one column , M+1 rows
  19. c = sym(c);                 % 符号类型矩阵
  20. b = ones(M+1,1);            % one column , M+1 rows
  21. delta=0.0001;   
  22. E1=9999; E0=0;
  23. t = 0;
  24. wep = 0.5 .* ones(M+1);     % M+1 square
  25. wes = 0.5 .* ones(M+1);

  26. while(1)
  27.     t = t + 1
  28.     % need revising
  29.     syms w
  30.     if(mode)
  31.         for i = 1:M+1           % length is odd
  32.             c(i) =  cos((i-1).* w) ;
  33.         end
  34.     else
  35.         for i = 1:M+1           % length is even
  36.              c = @(w) cos((i-1/2) .* w) ;
  37.         end
  38.     end
  39.         epw =  (1-c)' .* b;
  40.         esw =  b' .*  c;
  41.         wep = wep .* abs(epw);
  42.         wes = wes .* abs(esw);
  43.         
  44.     if(mode)
  45.         for n = 1:M+1           % length is odd
  46.             for m = n:M+1
  47. %                 Ppef = ( 1-cos((n-1)*w )).*( 1-cos((m-1)*w) );
  48. %                 Psef = (cos((n-1)*w ) .*cos((m-1)*w ));
  49. %                 Ep = Ppef(n) * wep(n) ;
  50. %                 Es = Psef(n) * wes(n) ;
  51.                 Ep =( 1-cos((n-1).*w )).*( 1-cos((m-1).*w) ) .* (wep(n,m) .* (abs((1-cos((n-1).* w)).*b(m))));   
  52.                 Epf = matlabFunction(Ep);
  53.                 Es = (cos((n-1).*w ) .*cos((m-1).*w )) .* wes(n,m) .* (abs((b(n).* c(m))));   
  54.                 Esf = matlabFunction(Es);
  55.                 P(n,m) = (quad(Epf,0,wp,0.000001) + quad(Esf,ws ,pi,0.000001 )) ;
  56.                 P(m,n) = P(n,m);
  57.             end
  58.         end
  59.     else
  60.         for n = 1:M+1           % length is even
  61.             for m = n:M+1
  62. %                 Ppef = @(w) ( 1-cos(((n-1)+1/2)*w )) .* ( 1-cos(((m-1)+1/2)*w) ) ;
  63. %                 Psef = @(w) cos(((n-1)+1/2)*w ) .* cos(((m-1)+1/2)*w );
  64.                 Ep =( 1-cos(((n-1)+1/2)*w )) .* ( 1-cos(((m-1)+1/2)*w) ) .* (wep(n,m) .* (abs((1-cos((n-1).* w)).*b(m))));
  65.                 Epf = matlabFunction(Ep);
  66.                 Es = (cos(((n-1)+1/2)*w ) .* cos(((m-1)+1/2)*w )) .* wes(n,m) .* (abs((b(n).* c(m))));   
  67.                 Esf = matlabFunction(Es);
  68.                 P(n,m) = (quad(Epf,0,wp,0.000001) + quad(Esf,ws ,pi,0.000001 )) ;
  69.                 P(m,n) = P(n,m);
  70.             end
  71.         end
  72.     end
  73.    
  74.     P;
  75.     [Vect,value] = eig(P,'nobalance');
  76.     b = Vect(:,1);
  77.    
  78.     % need revising
  79.     if(mode)
  80.         for i = 1:M             % length is odd
  81.             h(i) = b(M+2-i)/2;
  82.             h(M+M+2-i) = b(M+2-i)/2;
  83.         end
  84.         h(M+1) = b(1);
  85.     else
  86.         for i=1:M               % length is even
  87.             h(i) = b(M+1-i)/2;
  88.             h(M+M+1-i) = b(M+1-i)/2;
  89.         end
  90.     end

  91.     E0 = max(abs(P))
  92.    
  93.     if ( abs(E1-E0) < delta )
  94.         break;
  95.     else
  96.         E1 = E0;
  97.     end
  98.    
  99. end

  100. % Plot filter frequency response
  101. fft_size = 1024;
  102. q = abs( fft(h,fft_size));
  103. f_range = [0:1/fft_size:0.5];
  104. figure,plot(f_range, 20*log10(q(1:fft_size/2+1)/max(q)));
  105. title( 'FIR Low-pass Eigenfilter ')
  106. ylabel( 'Magnitude Response in(dB) ');
  107. xlabel( 'Normalized Frequency (f)');

  108. set(gca,'YTick',[ -80,-64,-48,-32,-16,0]);
  109. set(gca,'XTick',[0,fp,fs ,0.5]);grid;
  110. axis([0 .5 -80 5]);         % need revising

  111. figure,plot(f_range, abs(q(1:fft_size/2+1)/max(q)));
  112. title( 'FIR Low-pass Eigenfilter ');
  113. ylabel( 'Magnitude Response ');
  114. xlabel( 'Normalized Frequency (f)');
  115. axis([0 fp 0.9 1.01]);      % need revising
复制代码

原理1

原理1

原理2

原理2

最佳答案


20141303 发表于 2022-6-14 21:07:23
仅供参考,出现该问题的原因在于,在当次循环中Es = (cos((n-1).*w ) .*cos((m-1).*w )) .* wes(n,m) .* (abs((b(n).* c(m)))); 代码运行结果使得Es 为0,进而使得转化为函数句柄时变成了@()0.0,即函数里没有了未知变量,建议加个判断,即为零时对他积分为零,不为零时正常积分。
回复此楼

4 条回复


20141303 发表于 2022-6-14 21:07:23
仅供参考,出现该问题的原因在于,在当次循环中Es = (cos((n-1).*w ) .*cos((m-1).*w )) .* wes(n,m) .* (abs((b(n).* c(m)))); 代码运行结果使得Es 为0,进而使得转化为函数句柄时变成了@()0.0,即函数里没有了未知变量,建议加个判断,即为零时对他积分为零,不为零时正常积分。
回复此楼

Alison_Ni 发表于 2022-6-15 21:44:35
20141303 发表于 2022-6-14 21:07
仅供参考,出现该问题的原因在于,在当次循环中Es = (cos((n-1).*w ) .*cos((m-1).*w )) .* wes(n,m) .* (a ...

谢谢,这个问题解决了,能再请教你一个问题吗,原理2最下面的这个式子要怎么实现呢?一个函数矩阵的迭代,还涉及到绝对值,不知道怎么对句柄函数取绝对值,对句柄函数矩阵做运算。麻烦了
  1.         for n = 1:M+1           % length is odd
  2.             for m = n:M+1
  3.                 c(n) = cos((n-1)* w) ;
  4.                 epw(n,m) =  (1-c(n)) * b(m);
  5.                 esw(n,m) = b(n) * c(m);
  6.                 fepw = matlabFunction(abs(epw(n,m)));
  7.                 fesw = matlabFunction(abs(esw(n,m)));
  8.                 wep(n,m) = wep(n,m) .* fepw;
  9.                 wes(n,m) = wes(n,m) .* fesw;

  10.             end
  11.         end
复制代码

Alison_Ni 发表于 2022-6-15 21:45:46
Alison_Ni 发表于 2022-6-15 21:44
谢谢,这个问题解决了,能再请教你一个问题吗,原理2最下面的这个式子要怎么实现呢?一个函数矩阵的迭代 ...

未定义与 'function_handle' 类型的输入参数相对应的运算符 '.*'。
出错 EigenFilter_Equiripple (line 108)
                wep(n,m) = wep(n,m) .* fepw;

20141303 发表于 2022-6-16 09:20:00
仅供参考
改成
  1. wep(n,m) = wep(n,m) .* fepw(w);
复制代码

试试
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表