用下程序设计双相移偶对称任意频率点全相位陷波器,结果如图一所示
close all;clc;clear all;
N=16;
w=ones(1,N);
w1=hanning(N)';
win1=conv(w,w1);
win=win1/sum(win1);
H1=[1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1];
h1=ifft(H1);
h11=[h1(2:end) h1];
g1=h11.*win*N;
v1=exp(2*pi*j*(-N+1:N-1)*(0.5-0.3)/N);
g11=g1.*v1;
[h,f]=freqz(g11,1,N*100,'whole');
mag1=abs(h);
subplot(4,1,1),plot(f,mag1,'r');
set(gca,'XTick',0:pi/N*2:pi*2)
set(gca,'XTickLabel',{0:N})
title(' 相移0.3后偶对称单边1陷波器振幅特性 N=16 ');
ylim([0 2])
xlim([0 2*pi])
grid
H2=[1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1];
h2=ifft(H2);
h22=[h2(2:end) h2];
g2=h22.*win*N;
v2=exp(2*pi*j*(-N+1:N-1)*(0.5+0.3)/N);
g22=g2.*v2;
[h,f]=freqz(g22,1,N*100,'whole');
mag2=abs(h);
subplot(4,1,2),plot(f,mag2,'r');
set(gca,'XTick',0:pi/N*2:pi*2)
set(gca,'XTickLabel',{0:N})
title(' 相移-0.3后偶对称单边2陷波器振幅特性 N=16');
ylim([0 2])
xlim([0 2*pi])
grid
g=g11+g22;
%g(N)=g(N)-1;%先用这一语的句计算修正直流项系数'corr',得2.7226e-005)
g(N)=g(N)-1-2.7226e-005;%减去'corr'值,得最终g(N)值
[h,f]=freqz(g,1,N*100,'whole');
mag3=abs(h)';
subplot(4,1,3),plot(f,mag3,'r');
set(gca,'XTick',0:pi/N*2:pi*2)
set(gca,'XTickLabel',{0:N})
title('组合后相移0.3偶对称陷波器振幅特性 N=16');
ylim([0 2])
xlim([0 2*pi])
grid
mag4=10*log(mag3);
subplot(4,1,4),plot(f,mag4,'r');
set(gca,'XTick',0:pi/N*2:pi*2)
set(gca,'XTickLabel',{0:N})
title(' 校正后相移0.3偶对称陷波器对数振幅特性 N=16');
xlim([0 2*pi])
ylim([-400 0])
grid
disp('计算修正直流项系数')
corr=mag3((3.5-0.3)*100+1);
图一 双相移偶对称任意频率点的全相位陷波器设计
在<数字信号全相位谱分析和滤波技术>一书p.196给出了偶对称任意点陷波器系数设计公式(9-52),它和p194的传统对称任意点陷波器系数设计公式(9-46)是相同的,这样,全相位任意点陷波器系数可用下程序实现,它不必ifft,不必移位, 陷波器系数即乘卷积窗的陷波信号本身取样值.
clc;close all;clear;
N=16;
n=-N+1:N-1;
m=3.2;
f=ones(1,N); b=hanning(N)'; win=conv(f,b); win=win/max(win);
g=-2/N.*win.*cos(2*m*n*pi/N);
g(N)=(N-2)/N;%先用这一语的句计算修正直流项系数'corr',得2.7226e-005
g(N)=(N-2)/N-2.7226e-005;%减去'corr'值,得最终g(N)值
[G,w]=freqz(g,1,N*100,'whole');
mag=abs(G);
subplot(2,1,1); plot(w,mag,'r');
title('(a)N=16陷波点3.2的全相位陷波器振幅特性')
set(gca,'XTick',0:pi/N*2:2*pi)
set(gca,'XTickLabel',{0:N})
xlim([0 2*pi])
grid on;
legend('简化公式所得频率曲线')
subplot(2,1,2); plot(w,20*log10(mag) ,'r');
title('(b)N=16陷波点3.2的全相位陷波器对数振幅特性')
set(gca,'XTick',0:pi/N*2:2*pi)
set(gca,'XTickLabel',{0:N})
xlim([0 2*pi])
grid on;
disp('计算修正直流项系数')
corr=mag(m*100+1)
图二 用简化公式N=16 m=3.2全相位陷波器设计
这二个方法设计的N=16 m=3.2全相位陷波器系数都等於
g =[-0.00049653 -0.00074646 0.0052514 0.010651 -0.0069624 -0.034315 -0.014808 0.050564 0.06236 -0.028023 -0.10247 -0.034559 0.095876 0.099173 -0.038474 0.87891 -0.038474 0.099173 0.095876 -0.034559 -0.10247 -0.028023 0.06236 0.050564 -0.014808 -0.034315 -0.0069624 0.010651 0.0052514 -0.00074646 -0.00049653]