音頻信號調制，解調，加噪，去噪，濾波，matlab實現

題目分析

調制解調通信系統實現。

錄制一段語音信號，對其進行幅度調制解調，對比兩端語音信號的頻譜，並播放看其是否有偏差。在此基礎上，將調制后的語音信號加上一個噪聲，再濾波，觀察恢復語音信號的失真現象。

需要完成以下任務：

• 采集一段語音信號
• 做原信號時域頻域圖
• 對語音信號進行幅度調制
• 做出解調后的頻譜
• 調制解調前后做對比
• 調制后加噪聲
• 解調觀察對比

實驗流程

原信號分析

本次實驗從網絡下載WAV格式的音頻信號並將其長度裁剪為5s方便處理。

一般音樂和語音信號都是雙聲道信號，時域和頻譜圖會有兩個顏色，所以要取單列來分析，通過x1=x(:,1)語句來實現。

首先利用audioread函數讀取音頻文件獲得其頻率等信息。

然后利用fft函數求傅里葉變換然后利用abs函數求幅值畫幅度譜，再用plot畫出頻譜圖。為了方便觀察，可以利用fftshift對其進行左右互換，實驗預留了一個接口is_shift當其值為1時整個實驗的圖像都采用fftshift變換，反之不做，方便靈活選擇。

最終做出時域和頻域的兩個圖形。

以上畫了頻譜圖，可以更直觀的看到信號的特征，包括截止頻率和音頻范圍，音頻信號主要集中在低頻段。

幅度調制

這里統一使用了頻率歸一化，便於從原信號中讀取截止頻率和設置載波頻率。

取一定頻率的載波對信號進行AM調制，變量am_f表示載波信號的頻率，經過反復調參，認為取0.6最為合適。這里取了0.6pi為調制載波頻率，與原信號相乘實現AM調制，這里用點乘轉置矩陣實現。

從圖中可以看出原音樂信號的截止頻率為0.2pi(圖片對稱軸減一)，這里設置了0.6pi對信號進行AM調制，原信號的調制相當於頻譜搬移，左移一個右移一個，調制的目的是便於信號在信道中傳輸。當調制頻率較高時，聲音響度低，幾乎只能聽見茲茲的聲音，信號幾乎完全失真，當調制頻率較低時，聲音很尖銳，響度較大，能聽出調子，但也有茲茲的聲音。

解調

• 首先對調制后的信號進行解調：將調制后的信號與調制時相同的載波相乘實現解調，這里用點乘轉置矩陣實現。

• 然后用巴特沃斯濾波器對解調信號進行濾波：首先求巴特沃斯濾波器的頻率響應，其中用到了buttord求滿足性能指標的濾波器階數N和3dB截止頻率wc、用butter計算模擬濾波器的傳輸函數Ha(s)、用freqz求頻響。然后用filter實現濾波。
• 最后播放調制、解調和濾波后的聲音：通過聲音變化感受調制、解調和濾波。

首先得到了解調后的時域波形和頻譜圖，可以看出解調后的信號並沒有完全恢復原信號，會夾雜一點調制過程中的載波，通過濾波后信號頻譜有了很大改善。播放聲音發現：巴特沃斯濾波后聲音清晰，基本和原來的音樂差不多，但是音樂稍微低沉。巴特沃斯濾波器的特點是通頻帶的頻率響應曲線平滑。

加噪

• 信號加隨機白噪聲：首先求得隨機白噪聲，幅值為0.05，用randn語句產生噪聲，然后將噪聲和之前求得的原信號調制后的信號疊加。
• 濾掉噪聲：我使用了巴特沃斯濾波器來濾噪，其中用到buttord求滿足性能指標的濾波器階數N和3dB截止頻率wc、用butter求s域的頻率響應的參數、用bilinear函數即利用雙線性變換實現頻率響應s域到z域的變換，然后用filter求濾波后信號。

從圖中可以看到隨機噪聲的頻譜圖以及信號隨機噪聲后的時域波形和頻譜，通過播放感受到了噪聲對信號的影響；
濾波之后對噪聲的改善很大，噪聲變小，原聲音更加清晰，只是巴特沃斯濾波會把一部分原信號頻率濾掉，聲音會有點低沉。

附：實驗代碼

clear;clc;close all;
[x,fs] = audioread('water.wav');%讀取音樂信號

%輸出頻率
fprintf('fs1: %i \n',fs )
%音樂語音信號分聲道處理
x_l=x(:,1);
n1=length(x_l);
t1=(0:(n1-1))/fs;%length取數列長度即元素個數

is_shift=1; % 預留接口，下面用來判斷頻譜是否使用fftshift

%=======================原音頻時域/頻域圖=================================%
figure('NumberTitle', 'off', 'Name', '原音頻時域/頻域波形圖');
%畫音樂信號時域圖
subplot(2,1,1);plot(t1,x_l);
axis([0,5,-0.5,0.5]);xlabel('時間t');ylabel('幅度');title('音樂時域波形');


%畫音樂信號頻域圖
x_l_fft=fft(x_l,n1);
if is_shift==1 x_l_fft=fftshift(x_l_fft);end;
f1=0:fs/n1:fs*(n1-1)/n1;

subplot(2,1,2);plot(f1,abs(x_l_fft));
axis([0,fs,0,200]);xlabel('頻率f');ylabel('幅度');title('音樂信號頻譜');

%============================AM調制=======================================%

n=0:n1-1;%所有元素
w=2*f1/fs;%歸一化
am_f=0.6;%調制載波頻率

%調制
x_cos=cos(n*am_f*pi);%時余弦信號
x_am=x_l.*x_cos';%調制信號
x_cos_fft=fft(x_cos); 
x_am_fft=fft(x_am);
if is_shift==1 x_am_fft=fftshift(x_am_fft);end;

%畫出載波和調制信號波形
figure('NumberTitle', 'off', 'Name', 'AM調制圖');
subplot(2,2,1);plot(t1,x_cos);
axis([0,0.005,-1,1]);title('余弦信號時域波形');xlabel('時間t');ylabel('幅度');

subplot(2,2,2);plot(w,abs(x_cos_fft));
title('余弦信號頻譜');xlabel('w');ylabel('幅度');
subplot(2,2,3);plot(t1,x_am);
axis([0,5,-0.5,0.5]);title('調制信號時域波形');xlabel('時間t');ylabel('幅度');

subplot(2,2,4);plot(w,abs(x_am_fft));
title('調制信號頻譜');xlabel('w');ylabel('幅度');

%播放調制信號
%sound(x_l,fs);pause(5);%播放原音樂
%sound(x_am,fs);pause(5);%低頻調制播放

%============================解調========================================%

x_am_r=x_am.*x_cos'; %解調信號
x_am_r_fft=fft(x_am_r); %傅里葉變換
if is_shift==1 x_am_r_fft=fftshift(x_am_r_fft);end;

%畫解調信號時域波形和頻譜
figure('NumberTitle', 'off', 'Name', '解調時域/頻域波形圖');
subplot(2,1,1);plot(t1,x_am_r);
axis([0,5,-0.5,0.5]);xlabel('時間t');ylabel('幅度');title('解調信號時域波形');

subplot(2,1,2);plot(w,abs(x_am_r_fft));
xlabel('w');ylabel('幅度');title('解調信號頻譜');


%巴特沃斯低通濾波器過濾信號，wp<ws為低通濾波器
wp=am_f;%設置通帶頻率
ws=am_f*1.1; %設置阻帶頻率
rp=1;%設置通帶波紋系數
rs=50; %設置阻帶波紋系數    

[N,wc]=buttord(wp,ws,rp,rs); %求滿足性能指標的濾波器階數N和3dB截止頻率wc
[b,a]=butter(N,wc); %計算模擬濾波器的傳輸函數Ha(s)
[Hd,Ww]=freqz(b,a);%求頻響

figure('NumberTitle', 'off', 'Name', '解調后濾波圖');
subplot(3,1,1);plot(Ww/pi,abs(Hd));
axis([0,1,0,1.5]);xlabel('w/π');ylabel('幅度');title('巴特沃斯頻率響應曲線');

%求巴特沃斯濾波器濾波信號
x_am_r_f=filter(b,a,x_am_r);%濾波后的音樂信號
x_am_r_f=x_am_r_f*2;

N3=length(x_am_r_f);
t1=(0:(N3-1))/fs;%時間
f1=0:fs/N3:fs*(N3-1)/N3;%頻率
w1=2*f1/fs;%歸一化
x_am_r_f_fft=fft(x_am_r_f);%傅里葉變換
if is_shift==1 x_am_r_f_fft=fftshift(x_am_r_f_fft);end;

subplot(3,1,2);plot(t1,x_am_r_f);
axis([0,5,-0.5,0.5]);
xlabel('時間t');ylabel('幅度');title('巴特沃斯濾波信號時域波形');

subplot(3,1,3);plot(w1,abs(x_am_r_f_fft));
xlabel('w');ylabel('幅度');title('巴特沃斯濾波信號頻譜');

%sound(x_l,fs);pause(5);%原聲音
%sound(x_am_r_f,fs);pause(5);%解調並用巴特沃斯濾波器濾波后的聲音 


%============================調制后信號加噪===================================%

%加隨機白噪聲
x_noise=0.05*randn(n1,1);%噪聲占比5%
x_noise_fft=fft(x_noise);%噪聲傅里葉變換
if is_shift==1 x_noise_fft=fftshift(x_noise_fft);end;

%X_noise=x1+x_noise;%對調制后信號加隨機噪聲 
X_noise=x_am+x_noise;%對調制后信號加隨機噪聲 
X_noise_fft=fft(X_noise);%加噪信號傅里葉變換
if is_shift==1 X_noise_fft=fftshift(X_noise_fft);end;

figure('NumberTitle', 'off', 'Name', '加噪濾波頻譜圖');

subplot(3,2,1);plot(w,abs(x_noise_fft));
xlabel('w');ylabel('幅度');title('隨機白噪聲頻譜圖');

subplot(3,2,2);plot(t1,X_noise);
axis([0,5,-0.5,0.5]);
xlabel('時間t');ylabel('幅度');title('加隨機噪聲信號時域波形');

subplot(3,2,3);plot(w,abs(X_noise_fft));
xlabel('w');ylabel('幅度');title('加隨機噪聲信號頻譜圖');


%=======================含噪信號解調濾波==================================%
%先對含噪調制信號解調
X_noise_r=X_noise.*x_cos'; %解調信號
X_noise_r_fft=fft(X_noise_r); %傅里葉變換
if is_shift==1 X_noise_r_fft=fftshift(X_noise_r_fft);end;

%畫解調信號時域波形和頻譜
subplot(3,2,4);plot(w,abs(X_noise_r_fft));
xlabel('w');ylabel('幅度');title('含噪調制信號解調頻譜');



%再用巴特沃斯濾波器去噪
noise_w = 0.4;
wp=noise_w;ws=noise_w+0.06;Rp=1;Rs=15;
[n1,Wn1]=buttord(wp,ws,Rp,Rs,'s');%求低通濾波器的階數和截止頻率
[bb1,aa1]=butter(n1,Wn1,'s');%求s域的頻率響應的參數
[bb2,aa2]=bilinear(bb1,aa1,1);%利用雙線性變換實現頻率響應s域到z域的變換
%%subplot(3,2,2);plot(abs(bb2),aa2);
%%xlabel('w/π');ylabel('幅度');title('巴特沃斯頻率響應曲線');

X_noise_r_f=filter(bb2,aa2,X_noise_r);%求濾隨機噪聲后的信號
%畫濾噪后的圖
filter_noisy_fft=fft(X_noise_r_f);
if is_shift==1 filter_noisy_fft=fftshift(filter_noisy_fft);end;

subplot(3,2,5);plot(t1,X_noise_r_f);
axis([0,5,-0.5,0.5]);
xlabel('時間t');ylabel('幅度');title('濾隨機噪聲后信號時域波形');

subplot(3,2,6);plot(w,abs(filter_noisy_fft));
xlabel('w');ylabel('幅度');title('濾隨機噪聲后信號頻譜');

%播放
%sound(x_l,fs);pause(5);%原音樂
sound(X_noise,fs);pause(5);%加隨機噪聲
sound(X_noise_r_f,fs);pause(5);%濾隨機噪聲