社区导航

 

搜索
查看: 409|回复: 1

[资料分享] farrow结构滤波器matlab实现源码与仿真

[复制链接]

1570

TA的帖子

0

TA的资源

一粒金砂(中级)

Rank: 2

发表于 2019-3-6 21:09 | 显示全部楼层 |阅读模式
matlab源程序如下:
clc;close all;clear


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


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


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


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


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


for i=1:Nfil
    d=delay(i);
    if d==0
        d=d+0.0000001;   % avoid sin(0)/0;
    end
    xvec(i)=d;
    if (N/2-round(N/2))==0
        D=d+N/2;
    else
        D=d+N/2-0.5;
    end;
    cT=zeros(N+1,1);
    p1=cT;
    cT(1)=wp;
    if round(D)==D
        p1(1)=wp;
    else
        p1(1)=sin(D*wp*pi)/(D*pi);
    end;
    for k=1:N   % compute the elements of the Toeplitz matrix (vector)
        kD=k-D;
        cT(k+1) =sin(k*wp*pi)/(k*pi);
        p1(k+1) =sin(kD*wp*pi)/(kD*pi);
    end;
    Pz=toeplitz(cT);
    h=Pz\p1;
    hvec(i,=h';      % store designed coeff. vector
end


for k=1:N+1
    cc=polyfit(xvec,hvec(:,k),P);  % fit P:th-order polynomial to each coeff set
    C(:,k)=cc';
end
for j=1:Nfil
    d=delay(j);
    if d==0
        d=d+0.0000001;   % add 0.001 to avoid sin(0)/0;
    end
    h = C(P+1,;        % coeffs. via pol. approximation
    for n=1:P
        h=h+d^n*C(P+1-n,;
    end
    h=h/sum(h);           % scale response at zero freq. to unity
    H = freqz(h,1,w*pi);
    magresp(j, = abs(H);
    uwphase=-unwrap(angle(H));
    phasdel(j, = uwphase(2:Npt)./(w(2:Npt).*pi); % avoid divide by zero
end


% figure(1);
% for i=1:Nfil
%   plot(1:L,hvec(i,,'k'); hold on
% end;
% xlabel('Time In Samples'); ylabel('Impulse Response');
%
% figure(2);
% for i = 1:Nfil
%   plot(w,magresp(i,,'k');hold on
% end;
% xlabel('Normalized Frequency'); ylabel('Magnitude');
%
% figure(3);
% for i = 1:Nfil
%   plot(w(2:Npt),phasdel(i,,'k'); hold on
% end;
% xlabel('Normalized Frequency'); ylabel('Phase Delay');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs = 16;  %采样频率
T = 1 / Fs;
F = 3.2;      %信号频率
N_Max = F*2^18;   %采样点数
f = F/Fs;


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




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


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




Data_in_1_Delay = 0;
Data_in_2_Delay = 0;


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


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




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


for i = 0 : N_Max/2-1
    Data_in_Mismatch(2*i+1) = Data_in_1_Mismatch(i+1);
    Data_in_Mismatch(2*i+2) = Data_in_2_Mismatch(i+1);
end


figure(1)
FFT_Analysis(Data_in_Mismatch,Fs)


figure(2)
FFT_Analysis(Data_in_Delay(end-2^15+1-50:end-50),Fs)


回复

使用道具 举报

1

TA的帖子

0

TA的资源

一粒金砂(初级)

Rank: 1

发表于 2020-5-8 21:42 | 显示全部楼层
请问这一行代码p1(1)=sin(D*wp*pi)/(D*pi);是什么意思?

回复

使用道具 举报

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

关闭

站长推荐上一条 1/8 下一条

  • 论坛活动 E手掌握

    扫码关注
    EEWORLD 官方微信

  • EE福利  唾手可得

    扫码关注
    EE福利 唾手可得

Archiver|手机版|小黑屋|电子工程世界 ( 京ICP证 060456 )

GMT+8, 2020-5-25 13:57 , Processed in 0.117694 second(s), 21 queries , Gzip On, MemCache On.

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