频谱校正与滤波器设计matlab程序

频谱校正与滤波器设计matlab程序1、数字滤波器Matlab程序代码如下：wc=(400/1000)*pi;w1=boxcar(81);w2=hamming(81);w3=hanning(81);n=1:1:81;hd=sin(wc*(n-41))./(pi*(n-41));hd(41)=wc/pi;h1=hd.*w1';h2=hd.*w2';h3=hd.*w3';[mag1,rad]=freqz(h1);[mag2,rad]=freqz(h2);[mag3,rad]=freqz(h3);figure(1);plot(rad,20*log10(abs(mag1)));xlabel('NormalizedFrequency(rad)');ylabel('NormaliaedMagnitude(dB)');axis([0,3,-80,0]);title('矩形窗滤波器');gridon;figure(2);plot(rad,20*log10(abs(mag2)));xlabel('NormalizedFrequency(rad)');ylabel('NormaliaedMagnitude(dB)');axis([0,3,-80,0]);title('hamming窗滤波器');gridon;figure(3);plot(rad,20*log10(abs(mag3)));xlabel('NormalizedFrequency(rad)');ylabel('NormaliaedMagnitude(dB)');axis([0,3,-80,0]);title('hanning窗滤波器');gridon;输出结果如下：2、频谱校正（比值校正法）Matlab程序如下：functionZ=bizhi(x,fs,N,nx1,nx2,method)[nx,mx]=size(x);ifmx==1,x=x';end[nx,mx]=size(x);ifmxNx=[xzeros(1,N-mx)];elsex=x(1:N);endw=hann(N,'periodic');ifmethod==2xf=fft(x.*w');xf=xf(1:N/2)/N*4;WindowType=2;elsexf=fft(x);xf=xf(1:N/2)/N*2;WindowType=1;endddf=fs/N;n1=fix(nx1/ddf);n2=round(nx2/ddf);A=abs(xf);[Amax,index]=max(A(n1:n2));index=index+n1-1;phmax=angle(xf(index));if(WindowType==1)indsecL=A(index-1)A(index+1);df=indsecL.*A(index-1)./(Amax+A(index-1))-(1-indsecL).*A(index+1)./(Amax+A(index+1));Z(1)=(index-1-df)*ddf;Z(2)=Amax/sinc(df);Z(3)=(phmax+pi*df);endif(WindowType==2)indsecL=A(index-1)A(index+1);df=indsecL.*(2*A(index-1)-Amax)./(Amax+A(index-1))-(1-indsecL).*(2*A(index+1)-Amax)./(Amax+A(index+1));Z(1)=(index-1-df)*ddf;Z(2)=(1-df^2)*Amax/sinc(df);Z(3)=(phmax+pi*df);endZ(3)=mod(Z(3),2*pi);Z(3)=Z(3)-(Z(3)pi)*pi;设计例子如下：Fs=1000;n=0:1/Fs:1;xn=10+15*sin(10*2*pi*n)+randn(size(n));nfft=1024;Xn=fftshift(fft(xn,nfft));Xn(513)=Xn(513)/2;A=2*abs(Xn)/Fs;subplot(2,1,1),plot(xn)subplot(2,1,2)stem([-nfft/2:nfft/2-1]*Fs/nfft,A)xlabel('频率'),ylabel('幅值')P=Xn.*conj(Xn)/nfft;xlim([0500])Z=bizhi(xn,Fs,nfft,8,12,2)输出结果如下这样得到Z为10.000315.10511.5726表示信号频率为10.0003，幅值为15.1051，初始相角为1.5726。初始相角是以余弦信号为准，正弦信号正好差pi/2。