找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 5111|回复: 4
收起左侧

farrow结构滤波器matlab实现源码与仿真

[复制链接]
ID:294744 发表于 2018-3-20 21:38 | 显示全部楼层 |阅读模式
这是编写的farrow结构的滤波器的matlab仿真
0.png

matlab源程序如下:
  1. clc;close all;clear

  2. N = 29; % filter order, odd better
  3. L = N+1;             % filter length;
  4. Npt = 256;           % no. of frequency points for plots
  5. w = (0:1:Npt-1)/Npt; % frequenc scan (0,1)

  6. delay = [0 0.1 0.2 0.3 0.4 0.5];   % delay range x=0..0.5
  7. Nfil = length(delay); % number of filters

  8. h = zeros(1,L);      % impulse response vector
  9. hvec=zeros(Nfil,L);  % impulse response coefficient matrix
  10. magresp = zeros(Nfil,Npt);
  11. phasdel = zeros(Nfil,Npt-1);
  12. xvec=zeros(Nfil,1);     % fractional delay vector

  13. P = 3; % polynomial order for FARROW structure (ca. 1-5)
  14. C=zeros(P+1,N+1);      % polynomial coeff. matrix

  15. wp = 0.85; % normalized bandwidth (0-1.0)

  16. for i=1:Nfil
  17.     d=delay(i);
  18.     if d==0
  19.         d=d+0.0000001;   % avoid sin(0)/0;
  20.     end
  21.     xvec(i)=d;
  22.     if (N/2-round(N/2))==0
  23.         D=d+N/2;
  24.     else
  25.         D=d+N/2-0.5;
  26.     end;
  27.     cT=zeros(N+1,1);
  28.     p1=cT;
  29.     cT(1)=wp;
  30.     if round(D)==D
  31.         p1(1)=wp;
  32.     else
  33.         p1(1)=sin(D*wp*pi)/(D*pi);
  34.     end;
  35.     for k=1:N   % compute the elements of the Toeplitz matrix (vector)
  36.         kD=k-D;
  37.         cT(k+1) =sin(k*wp*pi)/(k*pi);
  38.         p1(k+1) =sin(kD*wp*pi)/(kD*pi);
  39.     end;
  40.     Pz=toeplitz(cT);
  41.     h=Pz\p1;
  42.     hvec(i,:)=h';      % store designed coeff. vector
  43. end

  44. for k=1:N+1
  45.     cc=polyfit(xvec,hvec(:,k),P);  % fit P:th-order polynomial to each coeff set
  46.     C(:,k)=cc';
  47. end
  48. for j=1:Nfil
  49.     d=delay(j);
  50.     if d==0
  51.         d=d+0.0000001;   % add 0.001 to avoid sin(0)/0;
  52.     end
  53.     h = C(P+1,:);        % coeffs. via pol. approximation
  54.     for n=1:P
  55.         h=h+d^n*C(P+1-n,:);
  56.     end
  57.     h=h/sum(h);           % scale response at zero freq. to unity
  58.     H = freqz(h,1,w*pi);
  59.     magresp(j,:) = abs(H);
  60.     uwphase=-unwrap(angle(H));
  61.     phasdel(j,:) = uwphase(2:Npt)./(w(2:Npt).*pi); % avoid divide by zero
  62. end

  63. % figure(1);
  64. % for i=1:Nfil
  65. %   plot(1:L,hvec(i,:),'k'); hold on
  66. % end;
  67. % xlabel('Time In Samples'); ylabel('Impulse Response');
  68. %
  69. % figure(2);
  70. % for i = 1:Nfil
  71. %   plot(w,magresp(i,:),'k');hold on
  72. % end;
  73. % xlabel('Normalized Frequency'); ylabel('Magnitude');
  74. %
  75. % figure(3);
  76. % for i = 1:Nfil
  77. %   plot(w(2:Npt),phasdel(i,:),'k'); hold on
  78. % end;
  79. % xlabel('Normalized Frequency'); ylabel('Phase Delay');

  80. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  81. Fs = 16;  %采样频率
  82. T = 1 / Fs;
  83. F = 3.2;      %信号频率
  84. N_Max = F*2^18;   %采样点数
  85. f = F/Fs;

  86. n1 = 0:2:(N_Max-1);
  87. n2 = 1:2:(N_Max-1);


  88. %*******加时间失配*******%
  89. Time_Mismatch = [0,-0.01];
  90. %Time_Mismatch = [0 0 0 0];

  91. Data_in_1_Mismatch = awgn( cos(2*pi*f*(n1+Time_Mismatch(1))),100 );
  92. Data_in_2_Mismatch = awgn( cos(2*pi*f*(n2+Time_Mismatch(2))),100 );


  93. Data_in_1_Delay = 0;
  94. Data_in_2_Delay = 0;

  95. % for delay_test = Time_Mismatch(2)-0.1:0.01:

  96. for i = 0:P
  97.     Data_in_1_Delay = Data_in_1_Delay + conv(C(P-i+1,:),Data_in_1_Mismatch)*(0.5*Time_Mismatch(1))^i;
  98.     Data_in_2_Delay = Data_in_2_Delay + conv(C(P-i+1,:),Data_in_2_Mismatch)*(0.5*Time_Mismatch(2))^i;
  99. end


  100. for i = 0 : N_Max/2 + N - 1
  101.     Data_in_Delay(2*i+1) = Data_in_1_Delay(i+1);
  102.     Data_in_Delay(2*i+2) = Data_in_2_Delay(i+1);
  103. end

  104. for i = 0 : N_Max/2-1
  105.     Data_in_Mismatch(2*i+1) = Data_in_1_Mismatch(i+1);
  106.     Data_in_Mismatch(2*i+2) = Data_in_2_Mismatch(i+1);
  107. end

  108. figure(1)
  109. FFT_Analysis(Data_in_Mismatch,Fs)

  110. figure(2)
  111. FFT_Analysis(Data_in_Delay(end-2^15+1-50:end-50),Fs)
复制代码

所有资料51hei提供下载:
Farrow.rar (3.08 KB, 下载次数: 62)
回复

使用道具 举报

ID:60531 发表于 2018-10-10 22:34 | 显示全部楼层
谢谢楼主分享
回复

使用道具 举报

ID:415698 发表于 2018-10-29 15:01 | 显示全部楼层
多谢楼主分享,希望有用
回复

使用道具 举报

ID:482797 发表于 2019-3-1 17:09 | 显示全部楼层
楼主,在吗? 最近在做个farrow滤波器,可以把你编程的思路的文档发我参考一下吗?看代码不是太清楚
回复

使用道具 举报

ID:482797 发表于 2019-3-1 17:09 | 显示全部楼层
楼主,在吗?
最近在做个farrow滤波器,可以把你编程的思路的文档发我参考一下吗?看代码不是太清楚
回复

使用道具 举报

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

本版积分规则

手机版|小黑屋|51黑电子论坛 |51黑电子论坛6群 QQ 管理员QQ:125739409;技术交流QQ群281945664

Powered by 单片机教程网

快速回复 返回顶部 返回列表