farrow結(jié)構(gòu)濾波器matlab實(shí)現(xiàn)源碼與仿真

ID:294744 · 發(fā)表于 2018-3-20 21:38

這是編寫的farrow結(jié)構(gòu)的濾波器的matlab仿真

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; %信號(hào)頻率
N_Max = F*2^18; %采樣點(diǎn)數(shù)
f = F/Fs;
n1 = 0:2:(N_Max-1);
n2 = 1:2:(N_Max-1);
%*******加時(shí)間失配*******%
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)