/
hanningFIR.m
158 lines (138 loc) · 3.79 KB
/
hanningFIR.m
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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
% 基于汉宁窗的FIR高通滤波器
clear
[x, Fs] = audioread("High_Frequency.wav");
sound(x,Fs);
pause(7);
% 原始信号初始分析
N=length(x); % 计算信号x的长度
t=0:1/Fs:(N-1)/Fs; % 计算时间范围,样本数除以采样频率
X=abs(fft(x));
X=X(1:N/2);
df=Fs/ N; % 计算频谱的谱线间隔
f=0:df:Fs/2-df; % 计算频谱频率范围
subplot(2,2,1);
plot(t,x);
title("原始音频信号时域图");
xlabel("时间(s)");
ylabel("幅度");
grid on;
subplot(2,2,2);
plot(f,X);
axis([900,1500,0,5000]);
title('原始音频信号频谱图');
xlabel('频率(hz)');
ylabel('幅度谱');
grid on;
% 加噪
fn1=930; % 单频噪声频率
fn2=950;
fn3=970;
fn4=1000;
noise=0.5*sin(fn1*2*pi*t)+0.5*sin(fn2*2*pi*t)+0.5*sin(fn3*2*pi*t)+0.5*sin(fn4*2*pi*t);
% 将噪声向量重复两次以形成 467097x2 的矩阵
noise_matrix = repmat(noise', 1, 2);
y=x+noise_matrix; % 加入一个多频噪声
sound(y,Fs);
pause(7);
% 加噪信号分析
Y=abs(fft(y));% 对加噪信号进行fft变换
Y=Y(1:N/2); % 截取前半部分
subplot(2,2,3);
plot(t,y);
title('加入四个单频噪声后的音频信号');
xlabel('时间(s)');
ylabel('幅度');
grid on;
subplot(2,2,4);
plot(f,Y);
axis([900,1500,0,5000]);
title('加入四个单频噪声后的音频信号频谱');
xlabel('频率(hz)');
ylabel('幅度谱');
grid on;
% 计算汉宁窗所需要的参数
fp=1100;fs=1020;%通带和阻带截止频率
Rp=1;As=43.9;%通带和阻带衰减
fcd=(fp+fs)/2;% 高通滤波器设计指标,截止频率
df=fp-fs; % 计算频率间隔
wc=fcd/Fs*2*pi;% 截止频率(弧度)
dw=df/Fs*2*pi;
wsd=fs/Fs*2*pi;
M=ceil(6.2*pi/dw)+1 % 计算汉宁窗设计该滤波器时需要的阶数
n=0:M-1; % 定义时间范围
w_ham=hann(M); % 产生M阶的汉宁窗
% M = 理想滤波器的长度
alpha = (M-1)/2;
n = 0:1:(M-1);
m = n - alpha + eps;
hd_bs=(sin(pi*m) - sin(wc*m))./(pi*m);
%%%%滤波器的形成
h_bs=w_ham'.*hd_bs; % 用窗口法计算实际滤波器脉冲响应
[H,w] = freqz(h_bs,1,1000,'whole');
H = (H(1:1:501))';
w = (w(1:1:501))';
mag = abs(H);
db = 20*log10((mag+eps)/max(mag));
pha = angle(H);
grd = grpdelay(h_bs,1,w);
% % % % 滤波器的性能指标
figure(2)
subplot(2,2,1);plot(w*Fs/(2*pi),db);
axis([0,1800,-120,20]);
title('幅度特性(db)');
xlabel('f(Hz)');ylabel('db');
grid on;
subplot(2,2,2);plot(w*Fs/(2*pi),mag);
axis([0,3900,-0.5,1.5]);
title('幅度特性曲线(Hz)');
xlabel('f(Hz)');ylabel('幅度');
grid on;
subplot(2,2,3);plot(w,pha);
axis([0,0.2,-3,3])
title('滤波器相频特性图');
xlabel('w/pi');ylabel('相位');
grid on;
subplot(2,2,4);plot(h_bs);
axis([250,600,-0.2,0.2]);
title('滤波器脉冲响应图');
xlabel('n');ylabel('h(n)');
grid on;
%%%%开始滤波
y_fil=filter(h_bs,1,y);% 用设计好的滤波器对y进行滤波
Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半
% % % 检验滤波
figure(3)
subplot(3,2,1);plot(t,x);
title('原音频信号时域');
xlabel('时间t');ylabel('幅度');
grid on;
subplot(3,2,2);plot(f,X);
axis([0,1500,0,3000]);
title('原音频信号幅度谱');
xlabel('频率f');ylabel('幅度');
grid on;
subplot(3,2,3);plot(t,y);
% axis([0,2.1,-0.05,0.1]);
title('加噪声后的音频信号时域');
xlabel('时间t');ylabel('幅度');
grid on;
subplot(3,2,4);plot(f,Y);
axis([0,1500,0,3000]);
title('加噪声后的音频信号幅度谱');
xlabel('频率f');ylabel('幅度');
grid on;
subplot(3,2,5);plot(t,y_fil);
% axis([0,2.1,-0.05,0.1]);
title('滤波后音频信号时域');
xlabel('时间t');ylabel('幅度');
grid on;
sound(y_fil,Fs)
subplot(3,2,6);plot(f,Y_fil);
axis([0,1500,0,3000]);
title('滤波后音频信号幅度谱');
xlabel('频率f');ylabel('幅度');
grid on;
figure(4)
zplane(h_bs,1);
xlabel('rez'); ylabel('jImz');
title('传输零极点');