image-20220210232819401

语音信号采集、时域、频域分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
%sound(y,fs);    %语音信号的播放
y = y(:,1); %取单声道
audiolen = length(y);
n = 0:audiolen-1;
dt = 1/fs;
t = n*dt;
disp(['语音信号长度为',num2str(max(t)),'s']);
figure('name','原语音信号分析');
%语音信号双声道时域波形
subplot(2,2,[1,2]);
plot(t,y);
xlabel('t/s');
ylabel('y');
title("语音信号时域波形图");
grid on;

Y = fft(y,audiolen);
N = length(Y);
df = fs/N;  %频域采样间隔
f = (0:N-1)*df; %频域采样点
%语音信号双声道幅频曲线
subplot(2,2,3);
plot(f,abs(Y));
xlabel('f/Hz');
ylabel('Y');
title("语音信号幅频曲线");
grid on;
%语音信号双声道相频曲线
subplot(2,2,4);
plot(f,angle(Y));
xlabel('f/Hz');
ylabel('angle(Y)/rad');
title("语音信号相频曲线");
grid on;

输出：

image-20211224201057154

时移

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y = y(:,1); %取单声道
dt = 1/fs;
%进行时移处理
y1 = [zeros(2*fs,1);y];
audiolen1 = length(y1);
t1 = (0:audiolen1-1)*dt;
audiowrite('时移2s音频.wav',y1,fs);
sound(y1,fs);   %播放时移后的信号
figure('name','经处理后的语音信号');
%语音信号时移2s后的双声道时域波形图
subplot(2,2,[1,2]);
plot(t1,y1);
xlabel('t/s');
ylabel('y');
title("经时移2s后的语音信号时域波形图");
grid on;

Y1 = fft(y1,audiolen1);
N1 = length(Y1);
df = fs/N1;  %频域采样间隔
f = (0:N1-1)*df; %频域采样点
%语音信号时移2s后的双声道幅频曲线
subplot(2,2,3);
plot(f,abs(Y1));
xlabel('f/Hz');
ylabel('Y');
title("语音信号幅频曲线");
grid on;

%语音信号时移2s后的双声道相频曲线
subplot(2,2,4);
plot(f,angle(Y1));
xlabel('f/Hz');
ylabel('angle(Y)/rad');
title("语音信号相频曲线");
grid on;

输出：

image-20211224201150402

单回声

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道
audiolen = length(y0);
n = 0:audiolen-1;
dt = 1/fs;
t = n*dt;

y1 = filter([1,zeros(1,30000),0.5],1,y0');
sound(5*y1,fs);

figure(1);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title("原信号时域波形图");

subplot(2,1,2);
plot(t,y1);
xlabel('t/s');
ylabel('y1');
title("单回声时域波形图");


df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
Y0 = fft(y0,audiolen);
audiolen1 = length(y1);
df1 = fs/audiolen1;  %频域采样间隔;
f1 = (0:audiolen1-1)*df1; %频域采样点
Y1 = fft(y1,audiolen1);

figure(2);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('Hz');
ylabel('Y0');
title("原信号幅频特性曲线");

subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('Hz');
ylabel('Y1');
title("单回声幅频特性曲线");

figure(3);
subplot(2,1,1);
plot(f,angle(Y0));
xlabel('Hz');
ylabel('Y0');
axis([0,100,-4,4]);
title("原信号相频特性曲线");

subplot(2,1,2);
plot(f1,angle(Y1));
xlabel('Hz');
ylabel('Y1');
axis([0,100,-4,4]);
title("单回声相频特性曲线");

image-20211224164324442

由时域波形图可知，时域波形上有较大差别，有的频率点原声幅度特别小但单回声由于延迟了所以在后面也还有很大幅度信号。

幅频曲线对比：

image-20211224165031331

单回响的幅度要比原信号大一些。

相频曲线对比：

image-20211224165554841

多重回声

以四重回声为例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道
N = 4;  %四重回声
audiolen = length(y0);
n = 0:audiolen-1;
dt = 1/fs;
t = n*dt;

y1 = filter(1,[1,zeros(1,80000/(N+1)),0.5],y0');
sound(5*y1,fs);

figure(1);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title("原信号时域波形图");

subplot(2,1,2);
plot(t,y1);
xlabel('t/s');
ylabel('y1');
title("四重回声时域波形图");


df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
Y0 = fft(y0,audiolen);
audiolen1 = length(y1);
df1 = fs/audiolen1;  %频域采样间隔;
f1 = (0:audiolen1-1)*df1; %频域采样点
Y1 = fft(y1,audiolen1);

figure(2);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('Hz');
ylabel('Y0');
title("原信号幅频特性曲线");

subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('Hz');
ylabel('Y1');
title("四重回声幅频特性曲线");

figure(3);
subplot(2,1,1);
plot(f,angle(Y0));
xlabel('Hz');
ylabel('Y0');
axis([0,100,-4,4]);
title("原信号相频特性曲线");

subplot(2,1,2);
plot(f1,angle(Y1));
xlabel('Hz');
ylabel('Y1');
axis([0,100,-4,4]);
title("四重回声相频特性曲线");

image-20211224170108704

image-20211224170120402

image-20211224170130888

无限回声

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道
audiolen = length(y0);
n = 0:audiolen-1;
dt = 1/fs;
t = n*dt;

a = 0.9;
Bz = [zeros(1,53),-a];
Az = [1,zeros(1,52),-a];
y1 = filter(Bz,Az,y0);
sound(5*y1,fs);

figure(1);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title("原信号时域波形图");

subplot(2,1,2);
plot(t,y1);
xlabel('t/s');
ylabel('y1');
title("无限回声时域波形图");


df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
Y0 = fft(y0,audiolen);
audiolen1 = length(y1);
df1 = fs/audiolen1;  %频域采样间隔;
f1 = (0:audiolen1-1)*df1; %频域采样点
Y1 = fft(y1,audiolen1);

figure(2);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('Hz');
ylabel('Y0');
title("原信号幅频特性曲线");

subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('Hz');
ylabel('Y1');
title("无限回声幅频特性曲线");

figure(3);
subplot(2,1,1);
plot(f,angle(Y0));
xlabel('Hz');
ylabel('Y0');
axis([0,100,-4,4]);
title("原信号相频特性曲线");

subplot(2,1,2);
plot(f1,angle(Y1));
xlabel('Hz');
ylabel('Y1');
axis([0,100,-4,4]);
title("无限回声相频特性曲线");

image-20211224170510637

image-20211224170527425

image-20211224170535179

全通混响

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道
audiolen = length(y0);
n = 0:audiolen-1;
dt = 1/fs;
t = n*dt;

a = 0.9;
Bz = [a,zeros(1,10),1];
Az = [1,zeros(1,10),a];
y1 = filter(Bz,Az,y0);
sound(5*y1,fs);

figure(1);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title("原信号时域波形图");

subplot(2,1,2);
plot(t,y1);
xlabel('t/s');
ylabel('y1');
title("全通混响时域波形图");


df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
Y0 = fft(y0,audiolen);
audiolen1 = length(y1);
df1 = fs/audiolen1;  %频域采样间隔;
f1 = (0:audiolen1-1)*df1; %频域采样点
Y1 = fft(y1,audiolen1);

figure(2);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('Hz');
ylabel('Y0');
title("原信号幅频特性曲线");

subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('Hz');
ylabel('Y1');
title("全通混响幅频特性曲线");

figure(3);
subplot(2,1,1);
plot(f,angle(Y0));
xlabel('Hz');
ylabel('Y0');
axis([0,100,-4,4]);
title("原信号相频特性曲线");

subplot(2,1,2);
plot(f1,angle(Y1));
xlabel('Hz');
ylabel('Y1');
axis([0,100,-4,4]);
title("全通混响相频特性曲线");

image-20211224170756974

image-20211224170806992

image-20211224170815881

通过全通滤波器的信号没有发生变化。

IIR滤波

低通

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道

%IIR滤波器参数
wp = 0.25*pi;
ws = 0.35*pi;
Fs = 5000;
Rp = 1; %通带最大衰减
As = 100;   %阻带最小衰减

T = 1/Fs;
wap = 2/T*tan(wp/2);
was = 2/T*tan(ws/2);

[N,wc] = buttord(wap,was,Rp,As,'s');
[Z,P,K] = buttap(N);
[bLP,aLP]=zp2tf(Z,P,K);
[bP,aP]=lp2lp(bLP,aLP,wc);
[b,a] = bilinear(bP,aP,Fs);
w = linspace(0,pi,100000);
H = freqz(b,a,w);
y1 = filter(b,a,y0);
Y0 = fft(y0);
Y1 = fft(y1);
sound(y1,fs);
figure(1);
plot(w/pi,abs(H));
xlabel('w/pi');
ylabel('H(jw)');
title('IIR低通滤波器频谱函数');

dt = 1/fs;
audiolen = length(y0);
n = 0:audiolen-1;
t = n*dt;
audiolen1 = length(y1);
n1 = 0:audiolen1-1;
t1 = n1*dt;

figure(2);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title('滤波前时域波形图');
subplot(2,1,2);
plot(t1,y1);
xlabel('n');
ylabel('y1');
title('滤波后时域波形图');

df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
df1 = fs/audiolen1;  %频域采样间隔
f1 = (0:audiolen1-1)*df1; %频域采样点
figure(3);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('f/Hz');
ylabel('Y0');
title('滤波前幅频特性');
subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('f/Hz');
ylabel('Y1');
title('滤波后幅频特性');

image-20211224190111103

image-20211224190423060

image-20211224190814433

可以看到高于875Hz左右的频率分量都被滤除了。

带通

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道

%带通滤波器参数
wpl = 0.2*pi;
wph = 0.75*pi;
wsl = 0.35*pi;
wsh = 0.85*pi;
Fs = 8000;
Rp = 1;
As = 100;
T = 1/Fs;
omepl = 2/T*tan(wpl/2);
omesl = 2/T*tan(wsl/2);
omesh = 2/T*tan(wsh/2);
omeph = 2/T*tan(wph/2);

%将模拟带通滤波器指标转换为归一化模拟低通滤波器的指标
Bp = omeph - omepl;
omep0 = sqrt(omepl*omeph);
omesl_ = ((omesl^2)-(omep0^2))/(omesl*Bp);
omesh_ = ((omesh^2)-(omep0^2))/(omesh*Bp);
omepl_ = ((omepl^2)-(omep0^2))/(omesh*Bp);
omeph_ = ((omeph^2)-(omep0^2))/(omeph*Bp);
if abs(omesh_)>abs(omesl_)  
    omes_ = abs(omesl_);
else
    omes_ = abs(omesh_);
end
omeap_ = 1;

[N,wc] = buttord(omeap_,omes_,Rp,As,'s');
[Z,P,K] = buttap(N);
[bLP,aLP] = zp2tf(Z,P,K);
[bBP,aBP] = lp2bp(bLP,aLP,omep0,Bp);
[b,a] = bilinear(bBP,aBP,Fs);
w = linspace(0,pi,100000);
H = freqz(b,a,w);
y1 = filter(b,a,y0);
Y0 = fft(y0);
Y1 = fft(y1);
sound(y1,fs);

figure(1);
plot(w,abs(H));
xlabel('w/pi');
ylabel('H(jw)');
title('IIR带通滤波器频谱函数');

dt = 1/fs;
audiolen = length(y0);
n = 0:audiolen-1;
t = n*dt;
audiolen1 = length(y1);
n1 = 0:audiolen1-1;
t1 = n1*dt;

figure(2);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title('滤波前时域波形图');
subplot(2,1,2);
plot(t1,y1);
xlabel('t/s');
ylabel('y1');
title('滤波后时域波形图');

df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
df1 = fs/audiolen1;  %频域采样间隔
f1 = (0:audiolen1-1)*df1; %频域采样点
figure(3);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('f/Hz');
ylabel('Y0');
title('滤波前幅频特性');
subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('f/Hz');
ylabel('Y1');
title('滤波后幅频特性');

image-20211224191037859

image-20211224190455197

image-20211224190502151

可以看出滤波后保留了1400Hz~3400Hz的频率分量，其他的都被滤除。

高通

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道

%IIR高通滤波器参数
wp = 0.7*pi;
ws = 0.9*pi;
Fs = 10000;
Rp = 1; %通带最衰减，As是阻带最小衰减
As = 100;
T = 1/Fs;
omep = (2/T)*tan(wp/2);
omes = (2/T)*tan(ws/2);
[N,wc] = buttord(omep,omes,Rp,As,'s');
[Z,P,K] = buttap(N);
[bLP,aLP] = zp2tf(Z,P,K);
G1 = tf(bLP,aLP);
[bHP,aHP] = lp2hp(bLP,aLP,wc);
[b,a] = bilinear(bHP,aHP,Fs);
w = linspace(0,pi,100000);
H = freqz(b,a,w);
y1 = filter(b,a,y0);
Y0 = fft(y0);
Y1 = fft(y1);
sound(y1,fs);


figure(1);
plot(w/pi,abs(H));
xlabel('w/pi');
ylabel('H(jw)');
title('IIR高通滤波器频谱函数');

dt = 1/fs;
audiolen = length(y0);
n = 0:audiolen-1;
t = n*dt;
audiolen1 = length(y1);
n1 = 0:audiolen1-1;
t1 = n1*dt;

figure(2);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title('滤波前时域波形图');
subplot(2,1,2);
plot(t1,y1);
xlabel('t/s');
ylabel('y1');
title('滤波后时域波形图');

df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
df1 = fs/audiolen1;  %频域采样间隔
f1 = (0:audiolen1-1)*df1; %频域采样点
figure(3);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('f/Hz');
ylabel('Y0');
title('滤波前幅频特性');
subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('f/Hz');
ylabel('Y1');
title('滤波后幅频特性');

image-20211224181629497

image-20211224191258554

image-20211224191305728

可看出低于截止频率的频率分量被滤除。

FIR滤波

高通

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道

%FIR高通滤波器参数
wp = 0.7*pi;
ws = 0.9*pi;
wc = 0.8*pi;
Fs = 8000;
Rp = 1;
As = 100;
wdel = ws-wp;
N = ceil(10*pi/wdel);
wn = (wp+ws)/2;

h = fir1(N,wn/pi,'high',kaiser(N+1));   %wn为归一化截止频率
y1 = conv(y0,h);
Y0 = fft(y0);
Y1 = fft(y1);
sound(5*y1,fs);


figure(1);
freqz(h,1,1024);
xlabel('w/pi');
ylabel('H(jw)');
title('FIR高通滤波器频谱函数');

dt = 1/fs;
audiolen = length(y0);
n = 0:audiolen-1;
t = n*dt;
audiolen1 = length(y1);
n1 = 0:audiolen1-1;
t1 = n1*dt;

figure(2);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title('滤波前时域波形图');
subplot(2,1,2);
plot(t1,y1);
xlabel('t/s');
ylabel('y1');
title('滤波后时域波形图');

df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
df1 = fs/audiolen1;  %频域采样间隔
f1 = (0:audiolen1-1)*df1; %频域采样点
figure(3);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('f/Hz');
ylabel('Y0');
title('滤波前幅频特性');
subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('f/Hz');
ylabel('Y1');
title('滤波后幅频特性');

image-20211224183153787

image-20211224191544088

image-20211224191553669

低通

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道

%FIR低通滤波器参数
wp = 0.25*pi;
ws = 0.3*pi;
Fs = 8000;
Rp = 1;
As = 100;
wdel = ws-wp;
N = ceil(10*pi/wdel);
wn = (wp+ws)/2;

h = fir1(N,wn/pi,'low',kaiser(N+1));   %wn为归一化截止频率
y1 = conv(y0,h);
Y0 = fft(y0);
Y1 = fft(y1);
sound(5*y1,fs);

figure(1);
freqz(h,1,1024);
xlabel('w/pi');
ylabel('H(jw)');
title('FIR低通滤波器频谱函数');
dt = 1/fs;
audiolen = length(y0);
n = 0:audiolen-1;
t = n*dt;
audiolen1 = length(y1);
n1 = 0:audiolen1-1;
t1 = n1*dt;

figure(2);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title('滤波前时域波形图');
subplot(2,1,2);
plot(t1,y1);
xlabel('t/s');
ylabel('y1');
title('滤波后时域波形图');

df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
df1 = fs/audiolen1;  %频域采样间隔
f1 = (0:audiolen1-1)*df1; %频域采样点
figure(3);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('f/Hz');
ylabel('Y0');
title('滤波前幅频特性');
subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('f/Hz');
ylabel('Y1');
title('滤波后幅频特性');

image-20211224183617166

image-20211224191831820

image-20211224191839215

带通

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道

%FIR带通通滤波器参数
wpl = 0.2*pi;
wsl = 0.35*pi;
wph = 0.75*pi;
wsh = 0.85*pi;
Fs = 8000;
Rp = 1;
As = 100;
wdel = abs(min((wsh-wph),(wpl-wsl)));
wp = (wpl+wsl)/2;
ws = (wph+wsh)/2;
wn = [wp,ws];
N = ceil(10*pi/wdel);

h = fir1(N,wn/pi,'bandpass',kaiser(N+1));   %wn为归一化截止频率
y1 = conv(y0,h);
Y0 = fft(y0);
Y1 = fft(y1);
sound(5*y1,fs);


figure(1);
freqz(h,1,1024);
xlabel('w/pi');
ylabel('H(jw)');
title('FIR带通滤波器频谱函数');
dt = 1/fs;
audiolen = length(y0);
n = 0:audiolen-1;
t = n*dt;
audiolen1 = length(y1);
n1 = 0:audiolen1-1;
t1 = n1*dt;

figure(2);
subplot(2,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title('滤波前时域波形图');
subplot(2,1,2);
plot(t1,y1);
xlabel('t/s');
ylabel('y1');
title('滤波后时域波形图');

df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
df1 = fs/audiolen1;  %频域采样间隔
f1 = (0:audiolen1-1)*df1; %频域采样点
figure(3);
subplot(2,1,1);
plot(f,abs(Y0));
xlabel('f/Hz');
ylabel('Y0');
title('滤波前幅频特性');
subplot(2,1,2);
plot(f1,abs(Y1));
xlabel('f/Hz');
ylabel('Y1');
title('滤波后幅频特性');

image-20211224201948485

image-20211224191942994

image-20211224191950298

思考题

image-20220210232949847

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
clear all;
close all;
clc;
[y,fs] = audioread('音频文件路径');   %语音信号采集
y0 = y(:,1); %取单声道
audiolen = length(y0);
noise = 0.05*randn(audiolen,1);
s = y0+noise;
%FIR低通滤波器参数
wp = 0.25*pi;
ws = 0.3*pi;
Fs = 8000;
Rp = 1;
As = 100;
wdel = ws-wp;
N = ceil(10*pi/wdel);
wn = (wp+ws)/2;

[b,a] = fir1(N,wn/pi,'low',kaiser(N+1));   %wn为归一化截止频率
y1 = filter(b,a,s);
Y0 = fft(y0);
Y1 = fft(y1);
S = fft(s);
sound(5*y1,fs);

dt = 1/fs;
audiolen = length(y0);
n = 0:audiolen-1;
t = n*dt;

figure(1);
subplot(3,1,1);
plot(t,y0);
xlabel('t/s');
ylabel('y0');
title('原始信号时域波形图');
subplot(3,1,2);
plot(t,s);
xlabel('t/s');
ylabel('y1');
title('加噪声后信号时域波形图');
subplot(3,1,3);
plot(t,y1);
xlabel('t/s');
ylabel('y1');
title('加噪声信号经FIR低通滤波后时域波形图');


df = fs/audiolen;  %频域采样间隔
f = (0:audiolen-1)*df; %频域采样点
figure(2);
subplot(3,1,1);
plot(f,abs(Y0));
xlabel('f/Hz');
ylabel('Y0');
title('原始信号幅频特性');
subplot(3,1,2);
plot(f,abs(S));
xlabel('f/Hz');
ylabel('S');
title('加噪声后信号幅频特性');
subplot(3,1,3);
plot(f,abs(Y1));
xlabel('f/Hz');
ylabel('Y1');
title('加噪声信号经FIR低通滤波后幅频特性');

image-20211224192304022

image-20211224192315970