数字信号处理实验

1.模拟信号的数字处理

clc;clear up;close all;

f0=13;

t=0:0.002:1;

x=cos(2*pi*f0*t);

subplot(2,2,1);

plot(t,x);

f1=10;

T1=1/f1;

nT1=0:T1:1;

x1=cos(2*pi*f0*nT1);

n=0:length(nT1)-1;

subplot(2,2,2);

stem(n,x1);

clc;clear up;close all;

f0=13;

t=0:0.002:1;

x=cos(2*pi*f0*t);

subplot(3,1,1);

plot(t,x);

f1=10;

T1=1/f1;

nT1=(0:T1:1)';

x1=cos(2*pi*f0*nT1);

n=0:length(nT1)-1;

subplot(3,1,2);

stem(n,x1);

t1 = linspace(-0.5,1.5,500)';

ya=sinc((1/T1)*t1(:,ones(size(nT1)))-(1/T1)*nT1(:,ones(size(t1)))')*x1;

subplot(3,1,3);

plot(t1,ya);

clc;

t=(0:0.002:1)';

xa=cos(80*pi*t)+sin(100*pi*t);

subplot(3,1,1);

plot(t,xa);

fs=200;

T=1/fs;

nT=(0:T:1)';

n=0:length(nT)-1;

x=cos(80*pi*nT)+sin(100*pi*nT);

subplot(3,1,2);

stem(n,x);

ya=sinc((1/T)*t(:,ones(size(nT))')-(1/T)*nT(:,ones(size(t)))')*x;

subplot(3,1,3)

plot(t,ya);

clc;

t=(0:0.002:1)';

xa=cos(80*pi*t)+sin(100*pi*t);

subplot(3,1,1);

plot(t,xa);

fs=150;

T=1/fs;

nT=(0:T:1)';

n=0:length(nT)-1;

x=cos(80*pi*nT)+sin(100*pi*nT);

subplot(3,1,2);

stem(n,x);

tic

% for i=1:size(t)

%     for j=1:size(nT)

%         a(j)=t(i)-nT(j);

%     end

%     ya(i)=sinc((1/T)*a)*x;

% end

% toc

% ya=zeros(1,100);

% n=1;

% while n<=fs+1

%     ya=ya+x(n)*sinc((t-nT(n))/T);

%     n=n+1;

% end

ya=sinc((1/T)*t(:,ones(size(nT))')-(1/T)*nT(:,ones(size(t)))')*x;

subplot(3,1,3)

plot(t,ya);

2离散信号的频域分析

clc;

x=[1,2,3,4];

w=linspace(-pi,pi,500);

[H,w]=freqz(x,1,w);

subplot(2,1,1);

plot(w/pi,abs(H));

subplot(2,1,2);

plot(w/pi,angle(H));

clc;

x=[1,2,3,4];

%DTFT

H=freqz(x,1,8,'whole');

figure;

subplot(2,1,1);

plot(0:7,abs(H));

subplot(2,1,2);

plot(0:7,angle(H));

%FFT

figure;

H1=fft(x,8);

subplot(2,1,1);

plot(0:7,abs(H1));

subplot(2,1,2);

plot(0:7,angle(H1));

clc;

fs=60;

T=1/fs;

nT=(0:T:1)';

n=0:length(nT)-1;

x=cos(80*pi*nT)+sin(100*pi*nT)+exp(1i*150*pi*nT);

y=fft(x);

plot(n,abs(y));

%若为10点DFT

% y=fft(x,10);

% plot(0:9,abs(y));

x1 = [1, 2, 3, 4, 2, 3, 4, 1];

x2 = [4, 3, 2, 1, 1, 2, 3, 2];

X1=fft(x1,8);

X2=fft(x2,8);

Y=ifft(X1.*X2);

C=cconv(x1,x2,8);

isequal(round(Y),C)

1. 离散时间系统的频域分析

close all;

clear up;

clc;

fs=2000;

T=1/fs;

nT=(0:T:1)';

n=0:length(nT)-1;

x1=[1 1 1 1];

x2=[1 2 3 4 4 3 2 1];

x3=[4 3 2 1 1 2 3 4];

x4=cos(pi/16*n);

x5=sin(pi/32*n);

x6=cos(1000*pi*nT)+cos(1200*pi*nT)+cos(2000*pi*nT);

x7=sin(100*pi*nT)+sin(200*pi*nT)+sin(300*pi*nT);

y1=fft(x1,4);

subplot(3,4,1);

stem(0:3,abs(y1));

subplot(3,4,2);

stem(0:3,angle(y1));

y2=fft(x2,8);

subplot(3,4,3);

stem(0:7,abs(y2));

subplot(3,4,4);

stem(0:7,angle(y2));

y3=fft(x3,8);

subplot(3,4,5);

stem(0:7,abs(y3));

subplot(3,4,6);

stem(0:7,angle(y3));

y4=fft(x4,8);

subplot(3,4,7);

stem(0:7,abs(y4));

subplot(3,4,8);

stem(0:7,angle(y4));

y5=fft(x5,8);

subplot(3,4,7);

stem(0:7,abs(y5));

subplot(3,4,8);

stem(0:7,angle(y5));

y6=fft(x5,8);

subplot(3,4,9);

stem(0:7,abs(y6));

subplot(3,4,10);

stem(0:7,angle(y6));

close all;

clear up;

clc;

fs=2000;

T=1/fs;

nT=(0:T:1)';

n=0:length(nT)-1;

x4=cos(pi/16*n);

x5=sin(pi/32*n);

x8=x4+x5;

y1=fft(x8,32);

subplot(3,4,1);

stem(0:31,abs(y1));

subplot(3,4,2);

stem(0:31,angle(y1));

close all;

clear up;

clc;

fs=2000;

T=1/fs;

nT=0:T:1;

n=0:length(nT)-1;

x1=[1 1 1 1];

x2=[1 2 3 4 4 3 2 1];

x3=[4 3 2 1 1 2 3 4];

x4=cos(pi/16*n);

x5=sin(pi/32*n);

x6=cos(1000*pi*nT)+cos(1200*pi*nT)+cos(2000*pi*nT);

x7=sin(100*pi*nT)+sin(200*pi*nT)+sin(300*pi*nT);

x8=x4+x5;

y6=fft(x6,555);

y7=fft(x7,555);

x = x6 + 1i * x7;  %构造信号

%% 计算16点x(n)的FFT

f = fft(x,555);

%% 利用对称性质求出x6(n)和x7(n)的傅里叶变换

N=555;

kk = conj(f); %取X的共轭

X6 = 0.5*(f+[kk(1),flip(kk(2:555))]);

X7 = -1i*0.5*(f-[kk(1),flip(kk(2:555))]);

figure;

plot((0:554), abs(X6),'r'); %取模，表示信号的频域幅值

figure;

plot((0:554), abs(y6),'r'); %取模，表示信号的频域幅值

4.FFT 的计算机实现

用三种不同的 DFT 程序计算 x(n) = R8 (n) 的傅里叶变换 X m ,并比较三种程序计算机运

行时间。

1. 用 for 语句的 M 函数文件 dft1.m，用循环变量逐点计算 X(k)；

dft1.m

 function[Am,pha]=dft1(x)

N=length(x);

w=exp(-j*2*pi/N);

for k=1:N

sum=0;

for n=1:N

sum=sum+x(n)*w^((k-1)*(n-1));

end

Am(k)=abs(sum);

pha(k)=angle(sum);

end

sl1.m

N=256;

x=[ones(1,8),zeros(1,N-8)];

t=cputime;

[Am1,pha1]=dft1(x);

t1=cputime-t

n=[0:(length(x)-1)];

w=(2*pi/length(x))*n;

figure(1)

subplot(2,1,1), plot(w,Am1,'b'); grid;

title('Magnitude part');

xlabel('frequency in radians');

ylabel('|X(exp(jw))|');

subplot(2,1,2), plot(w,pha1,'r'); grid;

title('Phase Part');

xlabel('frequency in radians');

ylabel('argX[exp(jw)]/radians');

1. 编写时域分解FFT 算法或频域分解FFT 算法的M函数文件dft2.m；

function [Am,pha]=dft2(x)

N=length(x); %输入数据的长度

M=log2(N); %蝶形运算的级数

indexs=zeros(1,N); %记录反序码值的数组

for n=0:N-1

    opposite=zeros(1,M); %记录反序码的数组

value=0; %记录反序码的值

    for m=1:M

        s=floor(n/2^(m-1)); %向下取整

        opposite(m)=mod(s,2); %计算余数

        value=value+opposite(m)*2^(M-m);

    end

    indexs(n+1)=value+1;

end

%%蝶形运算

%%W^J*(2^M-L)

%%M=log2(N)

%%J=0,1,2...(2^(L-1))

%%L=1,2 ,3...M

sz=zeros(1,N); %定义输出数组

for n=1:N

   sz(n)=x(indexs(n));

end

for L=1:M %按照 蝶形运行的级数依次计算

    d=2^(L-1); %两个点的距离

    for J=0:d-1

        p=J*(2^(M-L));

        W=exp(-i*2*pi*p/N); %转换因子公式，其中i代表虚部

        for K=J+1:2^L:N %找到每一个蝶形的左上点

            t1=sz(K);

            t2=sz(K+d);

            sz(K)=t1+t2*W; %蝶形运算公式

            sz(K+d)=t1-t2*W;

        end

    end

    Am=abs(sz);

    pha=angle(sz);

end

Xh.m

clc

close

clear

fs=2300;

N=4096;

n=0:N-1;

x=cos(1100*pi*n/fs)+cos(1300*pi*n/fs)+cos(2000*pi*n/fs);

t=cputime;

[Am,pha]=dft2(x);

t2=cputime-t

n=[0:(length(x)-1)];

w=(2*pi/length(x))*n;

figure(2)

subplot(211);

plot(w,Am,'b');grid;

title('Magnitude part');

xlabel('frequency in radians');

ylabel('|X(exp(jw))|');

subplot(212);

plot(w,pha,'r');grid;         %绘图

title('Phase Part');

xlabel('frequency in radians');

ylabel('argX[exp(jw)]/radians');

3、调用FFT库函数，直接计算X(k)；

Dft3.m

function[Am,pha]=dft3(x)

Xk=fft(x);

Am=abs(Xk);

pha=angle(Xk);

sl2.m

N=256;

x=[ones(1,8),zeros(1,N-8)];

t=cputime;

[Am3,pha3]=dft3(x);

t3=cputime-t

n=[0:(length(x)-1)];

w=(2*pi/length(x))*n;

figure(3)

subplot(2,1,1), plot(w,Am3,'b'); grid;

title('Magnitude part');

xlabel('frequency in radians');

ylabel('|X(exp(jw))|');

subplot(2,1,2), plot(w,pha3,'r'); grid;

title('Phase Part');

xlabel('frequency in radians');

ylabel('argX[exp(jw)]/radians');

设计IIR 数字滤波器

clc;clear;close all;

T=0.2*pi;Fs=1/T;

wpz=0.2;    

wsz=0.3;    

rp = 1;       

rs = 15;       

wp=2/T*tan(wpz*pi/2);ws=2/T*tan(wsz*pi/2);

[N, wc] = buttord(wp, ws, rp, rs, 's');

[B, A] = butter(N, wc, 'low','s');

[Bz, Az] = bilinear(B, A,Fs);  

[H, w] = freqz(Bz, Az,25,pi);

wp1 = 0.2;

ws1 = 0.3;

Rp1 = 2;

Rs1 = 30;

WP=2/T*tan(wp1*pi/2);WS=2/T*tan(ws1*pi/2);

[n, wc] = cheb1ord(WP,WS, Rp1, Rs1, 's');

[b, a] = cheby1(n, Rp1, wc,'low', 's');

[Bz1, Az1] = bilinear(b, a,Fs);  

[H1, w1] = freqz(Bz1, Az1,25,pi);

figure

subplot(2,2,1)

plot(w/pi,20*log10(abs(H)));

xlabel('\omega/pi');

ylabel('幅度/(dB)');

title('巴特沃斯低通IIR数字滤波器幅频响应特性曲线');

grid on;

subplot(2,2,2)

plot(w/pi,angle(H));

xlabel('\omega/pi');

ylabel('相位');

title('巴特沃斯低通IIR数字滤波器幅频响应特性曲线');

grid on;

subplot(2,2,3)

plot(w1/pi,20*log10(abs(H1)));

xlabel('\omega/pi');

ylabel('幅度/(dB)');

title('切比雪夫低通IIR数字滤波器幅频响应特性曲线');

grid on;

subplot(2,2,4)

plot(w1/pi,angle(H1));

xlabel('\omega/pi');

ylabel('相位');

title('切比雪夫低通IIR数字滤波器幅频响应特性曲线');

grid on;

clc;clear;close all;

T=0.2*pi;Fs=1/T;

wpz=0.2;    

wsz=0.3;    

• rp = 1;       

rs = 15;       

wp=2/T*tan(wpz*pi/2);ws=2/T*tan(wsz*pi/2);

[N, wc] = buttord(wp, ws, rp, rs, 's');

[B, A] = butter(N, wc, 'low','s');

[Bz, Az] = bilinear(B, A,Fs);  

[H, w] = freqz(Bz, Az,25,pi);

wp1 = 0.2;

ws1 = 0.3;

Rp1 = 2;

Rs1 = 30;

WP=2/T*tan(wp1*pi/2);WS=2/T*tan(ws1*pi/2);

[n, wc] = cheb1ord(WP,WS, Rp1, Rs1, 's');

[b, a] = cheby1(n, Rp1, wc,'low', 's');

[Bz1, Az1] = bilinear(b, a,Fs);  

[H1, w1] = freqz(Bz1, Az1,25,pi);

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,...

-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,...

-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

 b= filter(Bz,Az, x);

 q= filter(Bz1,Az1,x);

n = 1:length(x);

subplot(3,1,1);

plot(n, x);

title('原始信号');

xlabel('采样点');

ylabel('幅度');

subplot(3,1,2);

plot(n, b);

title('巴特沃斯滤波后的信号');

xlabel('采样点');

ylabel('幅度');

subplot(3,1,3);

plot(n, q);

title('切比雪夫滤波后的信号');

xlabel('采样点');

ylabel('幅度');

clc;clear;close all;

Fs=44100;

wpz=0.2;    

wsz=0.3;    

rp = 1;       

rs = 15;       

wp=2*Fs*tan(wpz*pi/2);ws=2*Fs*tan(wsz*pi/2);

[N, wc] = buttord(wp, ws, rp, rs, 's');

[B, A] = butter(N, wc, 'low','s');

[Bz, Az] = bilinear(B, A,Fs);  

[H, w] = freqz(Bz, Az,25,pi);

wp1 = 0.2;

ws1 = 0.3;

Rp1 = 2;

Rs1 = 30;

WP=2*Fs*tan(wp1*pi/2);WS=2*Fs*tan(ws1*pi/2);

[n, wc] = cheb1ord(WP,WS, Rp1, Rs1, 's');

[b, a] = cheby1(n, Rp1, wc,'low', 's');

[Bz1, Az1] = bilinear(b, a,Fs);  

[H1, w1] = freqz(Bz1, Az1,25,pi);

load('xas.mat');

 b= filter(Bz,Az, xas);

  q= filter(Bz1,Az1,xas);

n = 1:length(xas);

subplot(3,1,1);

plot(n, real(xas));

%sound( real(xas),Fs);

title('原始语言信号');

xlabel('采样点');

ylabel('幅度');

subplot(3,1,2);

plot(n, real(b));

sound( real(b),Fs);

title('巴特沃斯滤波后的语言信号');

xlabel('采样点');

ylabel('幅度');

subplot(3,1,3);

plot(n, real(q));

title('切比雪夫滤波后的语言信号');

xlabel('采样点');

ylabel('幅度');

%sound( real(q),Fs);

用窗函数法设计 FIR 数字滤波器

% 产生窗函数并画图

N = 64; % 窗长度

n = 0:N-1;

w_rect = ones(1,N);

w_hann = 0.5*(1-cos(2*pi*n/(N-1)));

w_hamming = 0.54-0.46*cos(2*pi*n/(N-1));

w_blackman = 0.42-0.5*cos(2*pi*n/(N-1))+0.08*cos(4*pi*n/(N-1));

figure;

plot(n,w_rect,'r',n,w_hann,'g',n,w_hamming,'b',n,w_blackman,'m');

title('Window Functions');

xlabel('Sample Index');

ylabel('Amplitude');

legend('Rectangular','Hanning','Hamming','Blackman');

% 在命令窗口输入以下命令即可生成并绘制矩形窗、汉宁窗、汉明窗和布莱克曼窗函数：

function [w] = window(type, N)

% type: 窗类型，可选 'rect', 'hanning', 'hamming', 'blackman'

% N: 窗长度

    switch type

        case 'rect'

            w = ones(N,1);

        case 'hanning'

            w = 0.5 - 0.5*cos(2*pi*(0:N-1)/(N-1));

        case 'hamming'

            w = 0.54 - 0.46*cos(2*pi*(0:N-1)/(N-1));

        case 'blackman'

            w = 0.42 - 0.5*cos(2*pi*(0:N-1)/(N-1)) + 0.08*cos(4*pi*(0:N-1)/(N-1));

        otherwise

            error('Unsupported window type!');

    end

    plot(w);

end

clc

close all

clear up

load('xas.mat');

% 滤波器参数

a = 0.4;

N1 = 15; % 窗长度为15

N2 = 33; % 窗长度为33

% 理想幅频特性

w = linspace(0, pi, 1000);

Hd = exp(-1j*w*a).*(w <= 0.4*pi);

% 矩形窗设计滤波器

h1_rect = fir1(N1-1, a/pi);

h2_rect = fir1(N2-1, a/pi);

H1_rect = freqz(h1_rect, 1, w);

H2_rect = freqz(h2_rect, 1, w);

[H1,w1]= freqz(h1_rect, 1, w);

[H2,w2]= freqz(h2_rect, 1, w);

% 汉宁窗设计滤波器

w_hann1 = hann(N1)';

w_hann2 = hann(N2)';

h1_hann = fir1(N1-1, a/pi, w_hann1);

h2_hann = fir1(N2-1, a/pi, w_hann2);

H1_hann = freqz(h1_hann, 1, w);

H2_hann = freqz(h2_hann, 1, w);

[H3,w3] = freqz(h1_hann, 1, w);

[H4,w4]= freqz(h2_hann, 1, w);

% 汉明窗设计滤波器

w_hamming1 = hamming(N1)';

w_hamming2 = hamming(N2)';

h1_hamming = fir1(N1-1, a/pi, w_hamming1);

h2_hamming = fir1(N2-1, a/pi, w_hamming2);

 H1_hamming = freqz(h1_hamming, 1, w);

H2_hamming = freqz(h2_hamming, 1, w);

[H5,w5]= freqz(h1_hamming, 1, w);

[H6,w6] = freqz(h2_hamming, 1, w);

%%

%语音信号恢复

Fs=44100;

load('xas.mat');

x = xas;

y = filter(h2_hamming, 1, x);

sound(real(y), Fs);

n = 1:length(xas);

figure(3);

subplot(3,1,1);

plot(n, real(xas));

subplot(3,1,2);

plot(n, real(y));

subplot(3,1,3);

plot(n, abs(fftshift(fft(y))));

%%

% 布莱克曼窗设计滤波器

w_blackman1 = blackman(N1)';

w_blackman2 = blackman(N2)';

h1_blackman = fir1(N1-1, a/pi, w_blackman1);

h2_blackman = fir1(N2-1, a/pi, w_blackman2);

H1_blackman = freqz(h1_blackman, 1, w);

H2_blackman = freqz(h2_blackman, 1, w);

[H7,w7] = freqz(h1_blackman, 1, w);

[H8,w8] = freqz(h2_blackman, 1, w);

% Kaiser窗设计滤波器

beta = 8.5;

w_kaiser1 = kaiser(N1, beta)';

w_kaiser2 = kaiser(N2, beta)';

h1_kaiser = fir1(N1-1, a/pi, w_kaiser1);

h2_kaiser = fir1(N2-1, a/pi, w_kaiser2);

H1_kaiser = freqz(h1_kaiser, 1, w);

H2_kaiser = freqz(h2_kaiser, 1, w);

[H9,w9] = freqz(h2_kaiser, 1, w);

[H10,w10] = freqz(h2_kaiser, 1, w);

% 绘图

figure(1);

subplot(2,2,1);

plot(w1/pi, 20*log10(abs(H1)), 'r', w3/pi, 20*log10(abs(H3)), 'g', w5/pi, 20*log10(abs(H5)), 'b', w7/pi, 20*log10(abs(H7)), 'm', w9/pi, 20*log10(abs(H9)), 'k');

title('损耗函数曲线(N=15)');

xlabel('\omega/pi');

ylabel('幅度/（dB）');

legend('Rectangular', 'Hanning', 'Hamming', 'Blackman', 'Kaiser');

subplot(2,2,2);

plot(w2/pi,20*log10(abs(H1)), 'r', w4/pi, 20*log10(abs(H4)), 'g', w6/pi, 20*log10(abs(H6)), 'b', w8/pi, 20*log10(abs(H8)), 'm', w10/pi, 20*log10(abs(H10)), 'k');

title('损耗函数曲线(N=33)');

xlabel('\omega/pi');

ylabel('幅度/（dB）');

legend('Rectangular', 'Hanning', 'Hamming', 'Blackman', 'Kaiser');

subplot(2,2,3);

plot(w1/pi, angle(H1), 'r', w3/pi, angle(H3), 'g', w5/pi, angle(H5), 'b', w7/pi, angle(H7), 'm', w9/pi, angle(H9), 'k');

title('相频特性曲线(N=15)');

xlabel('w/pi');

ylabel('Phase (rad)');

legend('Rectangular', 'Hanning', 'Hamming', 'Blackman', 'Kaiser');

subplot(2,2,4);

plot(w2/pi, angle(H2), 'r', w4/pi, angle(H4), 'g', w6/pi, angle(H6), 'b', w8/pi, angle(H8), 'm', w10/pi, angle(H10), 'k');

title('相频特性曲线(N=33)');

xlabel('w/pi');

ylabel('Phase (rad)');

legend('Rectangular', 'Hanning', 'Hamming', 'Blackman', 'Kaiser');

figure(2);

subplot(1,2,1);

plot(w/pi, abs(H1_rect), 'r', w/pi, abs(H1_hann), 'g', w/pi, abs(H1_hamming), 'b', w/pi, abs(H1_blackman), 'm', w/pi, abs(H1_kaiser), 'k');

title('损耗函数曲线(N=15)');

xlabel('\omega/pi');

ylabel('幅度/（dB）');

legend('Rectangular', 'Hanning', 'Hamming', 'Blackman', 'Kaiser');

subplot(1,2,2);

plot(w/pi, abs(H2_rect), 'r', w/pi, abs(H2_hann), 'g', w/pi, abs(H2_hamming), 'b', w/pi, abs(H2_blackman), 'm', w/pi, abs(H2_kaiser), 'k');

title('损耗函数曲线(N=33)');

xlabel('\omega/pi');

ylabel('幅度/（dB）');

legend('Rectangular', 'Hanning', 'Hamming', 'Blackman', 'Kaiser');

%%

%语音信号恢复

Fs=44100;

load('xas.mat');

x = xas;

y = filter(h2_hamming, 1, x);

sound(real(y), Fs);

n = 1:length(xas);

subplot(3,1,1);

plot(n, real(xas));

subplot(3,1,2);

plot(n, real(y));

subplot(3,1,3);

plot(n, abs(fftshift(fft(y))));

%%

x=[-4 -2 0 -4 -6 -4 -2 -4 -6 -6 -4 -4 -6 -6 -2 6 12 8 0 -16 -38 -60 -84 -90 66 -32 -4 -2 -4 8 12 12 10 6 6 6 4 0 0 0 0 0 -2 -4 0 0 0 -2 -2 0 0 -2 -2 -2 -2 0];

 b= filter(h2_hamming, 1,x);

n = 1:length(x);

figure(4)

subplot(2,1,1);

plot(n,x);

title('原始信号');

xlabel('采样点');

ylabel('幅度');

subplot(2,1,2);

plot(n, b);

title('汉宁窗设计滤波器滤波后的信号');

%语音信号恢复

Fs=44100;

 load('xas.mat');

 x = xas;

% y = filter(h2_hamming, 1, x);

% sound(real(y), Fs);

% n = 1:length(xas);

% figure(3);

% subplot(3,1,1);

% plot(n, real(xas));

% subplot(3,1,2);

% plot(n, real(y));

% subplot(3,1,3);

% plot(n, abs(fftshift(fft(y))));

c1=filter(h2_rect,1,xas);

d1=filter(h2_hann,1,xas);

e1= filter(h2_hamming, 1,xas);

f1=filter(h2_blackman,1,xas);

g1=filter(h2_kaiser,1,xas);

n = 1:length(x);

figure

subplot(3,4,1);

plot(n,real(xas));

subplot(3,4,2);

plot(n, abs(fftshift(fft(xas))));

title('原始信号');

xlabel('采样点');

ylabel('幅度');

subplot(3,4,3);

plot(n, real(c1));

subplot(3,4,4);

plot(n, abs(fftshift(fft(c1))));

title('Rectangular设计滤波器滤波后的信号');

subplot(3,4,5);

plot(n, real(d1));

subplot(3,4,6);

plot(n, abs(fftshift(fft(d1))));

title('Hanning设计滤波器滤波后的信号');

xlabel('采样点');

ylabel('幅度');

subplot(3,4,7);

plot(n, real(e1));

subplot(3,4,8);

plot(n, abs(fftshift(fft(e1))));

title('Hamming设计滤波器滤波后的信号');

subplot(3,4,9);

plot(n, real(f1));

subplot(3,4,10);

plot(n, abs(fftshift(fft(f1))));

title('Blackman设计滤波器滤波后的信号');

xlabel('采样点');

ylabel('幅度');

subplot(3,4,11);

plot(n, real(g2));

subplot(3,4,12);

plot(n, abs(fftshift(fft(g2))));

title('Kaiser设计滤波器滤波后的信号');

sound(real(g1), Fs);

fp=2*pi*1.5*10^3;fs=2*pi*3*10^3;rs=50;Fs=2*pi*1.5*10^4;

wp=2*pi*fp/Fs;

ws=2*pi*fs/Fs;

Bt=ws-wp;

alph=0.5842*(rs-21)^0.4+0.07886*(rs-21);

M=ceil((rs-8)/2.285/Bt);

wc=(wp+ws)/2/pi;

hn=fir1(M,wc,kaiser(M+1,alph));

%% 绘制幅频响应曲线

[H, w] = freqz(hn, 1, 512, Fs);

figure;

subplot(2,1,1);

plot(w, 20*log10(abs(H)));

title('幅频响应曲线');

xlabel('频率(rad/sec)');

ylabel('幅频(dB)');

%% 绘制相频响应曲线

subplot(2,1,2);

plot(w, angle(H));

title('相频响应曲线');

xlabel('频率(rad/sec)');

ylabel('相位(rad)');

wpl=0.55*pi;

wpu=0.7*pi;

wsl=0.45*pi;

wsu=0.8*pi;rs=40;

wc=[(wpl+wsl)/2,(wpu+wsu)/2];

Bt=wpl-wsl;

N=ceil(6.2*pi/Bt);

hn=fir1(N-1,wc/pi,hanning(N));

subplot(2,2,1);yn='h(n)';

stem(hn);

xlabel('n');

ylabel('h(n)');

subplot(2,2,2);A=1;

[H,w]=freqz(hn);

plot(w/pi,20*log10(abs(H)));

xlabel('\omega/pi');

ylabel('幅度/(dB)');

title('损耗函数曲线');