这里我们在Matlab中进行一个语音信号均匀量化的实验，关于量化的原理，这里不在详述，我们直接讲怎么进行语音量化，
Matlab版本：R2016a
那么我们的程序流程图如下：
图1 程序流程图
我们对程序流程图的解读如下：
读入语音信号这里我们是用的Matlab函数原型是[y,Fs] = audioread(filename)，其中filename是我们需要读入的文件名称，这里最好给读入文件的绝对地址，y,Fs分别输出的数据与采样频率。
画出读入语音信号的波形图，这里我们使用的函数原型是plot(X,Y)，这里我们将Y设置为读入的数据，X设置为Y的长度数组即可，那么我们画出的波形如下：
图2 读入语音时域图
找出语音中的最大值与最小值，这里我们选用的Matlab中的函数原型是[M,I] = max(___)，[M,I] = min(___)，这里输出的M代表的是最大值或者最小值，I代表的是最大值或者最小值在数组中的位置。
这里我们定义的量化个数为65536个即2的16次方。然后我们需要计算最小量化区间即(Max-Min)/M ，其中Max 和Min 分别为语音信号中的最大值与最小值，M是我们定义的量化个数，这里为65536，那么我们就可以计算得到相应的量化区间。
进行量化，这里我们将当前语音信号值减去最小值，然后在除以最小量化间隔，这样就可以得到相应的量化值了，然后我们再进行解量化，即我们用相应的量化值乘以量化间隔，然后再加上最小值，即得源信号值，这里我们将解量化后的值与源信号进行对比，即取两个信号差值的绝对值，得到的图像如下：
图3 差值图
我们会发现，相应的差值基本一样，即这个插值为量化误差。然后我们要进行相应的调制，这里我们首先需要将量化的十进制转化为并行的16进制，这里我们先做除法，然后取余数即可，这里用到的Matlab函数原型为：b = mod(a,m)，其中a为被除数，m为除数，b为余数。将其转换为并行的二进制之后，我们需要将其变换为串行的二进制数，这里用到的Matlab函数原型为B = reshape(A,sz1,...,szN)，这里A为输入的数组，后面的数字即转换为多少行与列，这里我们将其转换为1列即可。
接下来我们要分别进行ASK调制与PSK调制，这里我们分别进行相应的讲解：
2ASK的程序流程图如下：
图4 2ASK流程图
添加高斯白噪声，这里我们选用的是awgn（x,snr,’measured’），这里x为需要添加噪声的序列，snr为需要的信噪比，’measured’代表的是计算输入序列的平均功率之后，按照信噪比进行添加；
抽样判决，这里抽样判决的最佳门限为b*=a/2 这里a代表的是输入的幅值大小，所以这里a设置为1，最佳噪声门限为0.5，当添加噪声的序列>0.5时我们判决为1，小于等于0.5时，我们判决为0；
进行误码率的计算，这里我们选用的是Matlab中的biterr函数，[number,ratio] = biterr(x,y)，进行比较的两个数组分别为x,y，经过比较后，numbers输出的是错误个数，ratio输出的是错误率，这样我们就得到了相应信噪比下的误码率，然后我们再进行多次循环，这样就能得到相应的误码率与信噪比的关系。
然后我们将不同信噪比下的波形时域图画出来，我们就可以观察到不同信噪比下的时域波形图：
图5 不同信噪比下ASK的解调、解量化后的时域波形图
我们会发现随着信噪比的增加，进行调制后解调的信号跟原信号的差距用肉眼已经观察不出来了，这里我们将其做差，求其差信号，其图如下：
图6 不同信噪比下ASK的解调、解量化后的时域波形与原波形的差值图
这里我我们可以很明显的看出，在15dB时其差值开始加快缩小，然后在20dB时，其误差基本可以忽略不计。
然后我们再进行PSK调制，PSK的程序流程图如下：
图7 2PSK程序流程图
进行2PSK的调制，这里我们需要将产生的随机序列进行简单的调制，当产生的随机序列为1时，我们将其调制为-1，当其产生的随机序列为0时，我们将其调制为1，这样就实现了相位调制。
添加高斯白噪声，这里我们是通过计算整个序列的平均功率，然后在平均功率的基础上进行噪声的添加。
抽样判决，这里抽样判决的最佳门限为b*=0 当添加噪声的序列>0时我们判决为0，小于等于0时，我们判决为1；
进行误码率的计算，这里我们选用的是Matlab中的biterr函数，[number,ratio] = biterr(x,y)，进行比较的两个数组分别为x,y，经过比较后，numbers输出的是错误个数，ratio输出的是错误率，这样我们就得到了相应信噪比下的误码率，然后我们再进行多次循环，这样就能得到相应的误码率与信噪比的关系。
然后我们将不同信噪比下的波形时域图画出来，我们就可以观察到不同信噪比下的时域波形图：
图8 不同信噪比下PSK的解调、解量化后的时域波形图
我们会发现随着信噪比的增加，进行调制后解调的信号跟原信号的差距用肉眼已经观察不出来了，这里我们将其做差，求其差信号，其图如下：
图9 不同信噪比下PSK的解调、解量化后的时域波形与原波形的差值图
我们发现当信噪比为15dB时其差值就已经趋近于一个定值了。说明信噪比已经基本不影响其误码率了，影响的另一个因素为量化误差。接下来我们画出了其实际的误码率曲线图，如下图：
图10 ASK、PSK调制误码率曲线图
从图中可以看出ASK的误码率随着信噪比的增加在下降，但是PSK的误码率在10dB之后却没有了，即其误码率为0，这说明PSK调制、解调后在该语音信号下没有发生差错，这从侧面也证实了在相同的信噪比下PSK的抗噪性能优于ASK的抗噪性能。
那么我们的Matalab的实现代码如下：
%语音信号的量化
%读取语音信号
[source,Fs] = audioread('E:\Matlab_Code\Quantification\voicefile.wav');
source_Quan_16bit = zeros(length(source),16);%定义16位数组
jungle_dec = zeros(length(source),1);
psk = zeros(16*length(source),1);%PSK调制数组
ration_ASK = zeros(5,1); %定义ASK调制误码率数组
ration_PSK = zeros(5,1); %定义PSK调制误码率数组
%得到最小值与最大值
[maxdata,I_Max] = max(source);
[mindata,I_Min] = min(source);
%设定量化区间
M = 65536;
v = (maxdata-mindata)/M;
%进行量化
source_Quan = uint16((source-mindata)./v);
%进行解量化：
source_DeQuan = double(source_Quan).* v + mindata;
%将数据写出来
audiowrite('E:\Matlab_Code\Quantification\my1.wav',source_DeQuan,Fs);
%规定时间区间
n=1:1:length(source);
%画出读入的信号
figure(1)
plot(n,source);
xlabel('时间t')
ylabel('幅值')
title('读入语音信号的时域图')
%对量化之后与源信号做差
source_difference = abs(source-double(source_DeQuan));
figure(2)
semilogy(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('原始信号与解量化之后的差值')
%进行并转串16位
for k=1:length(source)
for z = 1:16
if(source_Quan(k) == 0)
source_Quan_16bit(k,z) = 0;
else
source_Quan_16bit(k,z) = mod(uint16(double(source_Quan(k))/2^(16-z)-0.5),2);
end
end
end
%将16bit转换为1bit
source_Quan_1bit = reshape(source_Quan_16bit,16*length(source),1);
%进行ASK调制
EbN0_dB = 0:5:20;
for i=1:length(EbN0_dB)
%添加高斯噪声
Bit_Noise = awgn(source_Quan_1bit,EbN0_dB(i),'measured');
%判决
jungle_1bit = Bit_Noise > 0.5;
%计算误码率
[number_ASK,ration_ASK(i)] = biterr(source_Quan_1bit,jungle_1bit);
%将数据进从1位变为16位
jungle_16bit = reshape(jungle_1bit,length(source),16);
%将16位2进制数变为10进制数
for k=1:length(source)
b = 0;
for j=1:16
b = b+jungle_16bit(k,j)*2^(16-j);
end
jungle_dec(k) = b;
end
%进行解量化,5个不同的信噪比下
if( i== 1)
jungle_DeQuan_1 = jungle_dec.* v + mindata;
end
if( i== 2)
jungle_DeQuan_2 = jungle_dec.* v + mindata;
end
if( i== 3)
jungle_DeQuan_3 = jungle_dec.* v + mindata;
end
if( i== 4)
jungle_DeQuan_4 = jungle_dec.* v + mindata;
end
if( i== 5)
jungle_DeQuan_5 = jungle_dec.* v + mindata;
end
end
%画图，画出不同信噪比下的波形图
figure(3)
subplot(2,3,1)
plot(n,source_DeQuan);
xlabel('time');
ylabel('幅值');
title('原始解量化信号')
subplot(2,3,2)
plot(n,jungle_DeQuan_1);
xlabel('time');
ylabel('幅值');
title('0dB信噪比下的信号图(ASK)')
subplot(2,3,3)
plot(n,jungle_DeQuan_2);
xlabel('time');
ylabel('幅值');
title('5dB信噪比下的信号图(ASK)')
subplot(2,3,4)
plot(n,jungle_DeQuan_3);
xlabel('time');
ylabel('幅值');
title('10dB信噪比下的信号图(ASK)')
subplot(2,3,5)
plot(n,jungle_DeQuan_4);
xlabel('time');
ylabel('幅值');
title('15dB信噪比下的信号图(ASK)')
subplot(2,3,6)
plot(n,jungle_DeQuan_5);
xlabel('time');
ylabel('幅值');
title('20dB信噪比下的信号图(ASK)')
%画出解调后信号与原信号的差值
figure(4)
%0dB
subplot(2,3,1)
source_difference = abs(source-double(jungle_DeQuan_1));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('ASK解调信号与原始信号的差值（0dB）')
%5dB
subplot(2,3,2)
source_difference = abs(source-double(jungle_DeQuan_2));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('ASK解调信号与原始信号的差值（5dB）')
%10dB
subplot(2,3,3)
source_difference = abs(source-double(jungle_DeQuan_3));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('ASK解调信号与原始信号的差值（10dB）')
%15dB
subplot(2,3,4)
source_difference = abs(source-double(jungle_DeQuan_4));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('ASK解调信号与原始信号的差值（15dB）')
%20dB
subplot(2,3,5)
source_difference = abs(source-double(jungle_DeQuan_5));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('ASK解调信号与原始信号的差值（20dB）')
%将数据写出来
audiowrite('E:\Matlab_Code\Quantification\my2.wav',jungle_DeQuan_5,Fs);
%进行PSK的调制
for k=1:length(source_Quan_1bit)
if(source_Quan_1bit(k) == 1)
psk(k) = -1;
elseif (source_Quan_1bit(k) == 0)
psk(k) = 1;
end
end
for i=1:length(EbN0_dB)
%添加噪声
Bit_Noise = awgn(psk,EbN0_dB(i),'measured');
%进行判决
jungle_1bit = Bit_Noise < 0;
%计算误码率
[number_PSK,ration_PSK(i)] = biterr(source_Quan_1bit,jungle_1bit);
%将数据进从1位变为16位
jungle_16bit = reshape(jungle_1bit,length(source),16);
%将16位2进制数变为10进制数
for k=1:length(source)
b = 0;
for j=1:16
b = b+jungle_16bit(k,j)*2^(16-j);
end
jungle_dec(k) = b;
end
%进行解量化,5个不同的信噪比下
if( i== 1)
jungle_DeQuan_1 = jungle_dec.* v + mindata;
end
if( i== 2)
jungle_DeQuan_2 = jungle_dec.* v + mindata;
end
if( i== 3)
jungle_DeQuan_3 = jungle_dec.* v + mindata;
end
if( i== 4)
jungle_DeQuan_4 = jungle_dec.* v + mindata;
end
if( i== 5)
jungle_DeQuan_5 = jungle_dec.* v + mindata;
end
end
%画图，画出不同信噪比下的波形图
figure(5)
subplot(2,3,1)
plot(n,source_DeQuan);
xlabel('time');
ylabel('幅值');
title('原始解量化信号')
subplot(2,3,2)
plot(n,jungle_DeQuan_1);
xlabel('time');
ylabel('幅值');
title('0dB信噪比下的信号图(PSK)')
subplot(2,3,3)
plot(n,jungle_DeQuan_2);
xlabel('time');
ylabel('幅值');
title('5dB信噪比下的信号图(PSK)')
subplot(2,3,4)
plot(n,jungle_DeQuan_3);
xlabel('time');
ylabel('幅值');
title('10dB信噪比下的信号图(PSK)')
subplot(2,3,5)
plot(n,jungle_DeQuan_4);
xlabel('time');
ylabel('幅值');
title('15dB信噪比下的信号图(PSK)')
subplot(2,3,6)
plot(n,jungle_DeQuan_5);
xlabel('time');
ylabel('幅值');
title('20dB信噪比下的信号图(PSK)')
%画出不同信噪比下解调信号与源信号的差值
figure(6)
%0dB
subplot(2,3,1)
source_difference = abs(source-double(jungle_DeQuan_1));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('PSK解调信号与原始信号的差值（0dB）')
%5dB
subplot(2,3,2)
source_difference = abs(source-double(jungle_DeQuan_2));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('PSK解调信号与原始信号的差值（5dB）')
%10dB
subplot(2,3,3)
source_difference = abs(source-double(jungle_DeQuan_3));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('PSK解调信号与原始信号的差值（10dB）')
%15dB
subplot(2,3,4)
source_difference = abs(source-double(jungle_DeQuan_4));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('PSK解调信号与原始信号的差值（15dB）')
%20dB
subplot(2,3,5)
source_difference = abs(source-double(jungle_DeQuan_5));
plot(n,source_difference)
xlabel('时间t')
ylabel('幅值差')
title('PSK解调信号与原始信号的差值（20dB）')
%将数据写出来
audiowrite('E:\Matlab_Code\Quantification\my3.wav',jungle_DeQuan_5,Fs);
%进行性能分析
figure(7)
semilogy(EbN0_dB,ration_ASK,'-or',EbN0_dB,ration_PSK,'-*b')
grid on
ylabel('Pe')
xlabel('r/dB')
legend('ASK误码率曲线图','PSK误码率曲线图');
完。
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
