数字代码

%信号的读入与频谱观察
name=input('Please input the name of the .wav file:','s');
[y,Fs,bits]=wavread(name);%读出信号，采样率和采样位数。
sigLength=length(y);
Y=fft(y,sigLength);
f=Fs*(1:sigLength)/sigLength;
absY=abs(Y);
figure(1)
subplot(3,2,1);plot(f,absY);xlabel('Frequency(Hz)');grid on
t=(0:sigLength-1)/Fs;
subplot(3,2,2);plot(t,y);xlabel('Time(s)');grid on
y1=y;
%第一级操作：
%内插1：
I=2;
ycz1=zeros(1,I*length(y1));
ycz1(1:I:length(ycz1))=y1;
yf1=fft(ycz1,length(ycz1));
f=(1:length(ycz1))/length(ycz1)*Fs*2;
figure;plot(f,abs(yf1),'k')
%滤波器设计
fs=8000;fp=7000;Fs=100000;Rp=1;Rs=90;
ws=2*fs/Fs;wp=2*fp/Fs;
[N,Wpo]=ellipord(wp,ws,Rp,Rs);
[b,a]=ellip(N,Rp,Rs,Wpo);
w=0:0.005*pi:pi;
[h,w]=freqz(b,a,w);
fy=fft(h,length(h));
figure;
plot(w/pi,20*log(abs(fy)));
%信号滤波
ycz12=filter(b,a,ycz1);
yf12=fft(ycz12,length(ycz12));
f=(1:length(ycz12))/length(ycz12)*Fs*2;
figure;plot(f,abs(yf12),'k')
%信号抽取；
D=5;
ycq1=ycz12(1:D:length(ycz12));
yf=fft(ycq1,length(ycq1));
f=(1:length(ycq1))/length(ycq1)*Fs*2/5;
figure;plot(f,abs(yf),'k');




%第二次内插
I=4;
ycz2=zeros(1,I*length(ycq1));
ycz2(1:I:length(ycz2))=ycq1;
yf2=fft(ycz2,length(ycz2));
f=(1:length(ycz2))/length(ycz2)*Fs*8/5;
figure;plot(f,abs(yf2),'k')
&滤波设计
fs=8000;fp=7000;Fs=80000;Rp=1;Rs=90;
ws=2*fs/Fs;wp=2*fp/Fs;
[N,Wpo]=ellipord(wp,ws,Rp,Rs);
[b,a]=ellip(N,Rp,Rs,Wpo);
w=0:0.005*pi:pi;
[h,w]=freqz(b,a,w);
fy=fft(h,length(h));
figure;
plot(w/pi,20*log(abs(fy)));
%信号滤波
ycz22=filter(b,a,ycz2);
yf22=fft(ycz22,length(ycz22));
f=(1:length(ycz22))/length(ycz22)*Fs*8/5;
figure;plot(f,abs(yf22),'k')
%信号抽取
D=5;
ycq2=ycz22(1:D:length(ycz22));
yf=fft(ycq2,length(ycq2));
f=(1:length(ycq2)/length(ycq2)*16000;
figure;plot(f,abs(yf),'k');