Matlab实现信号分析与处理

单位样本序列生成

$$\delta [n]=\left\{\begin{matrix} 1,& n=0\\ 0,&n\ne 0 \end{matrix}\right.$$

% 横坐标 -10 to 20
n = -10:20;
% 纵坐标，与横坐标相对应
x = [zeros(1,10) 1 zeros(1,20)];
% 生成图像
stem(n,x);
xlabel('n');ylabel('\delta[n]'); 

时域线性卷积

$$y[n]=x[n]*h[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]$$

h = [3 2 1 -2 1 0 -4 0 3];  % 脉冲响应
x = [1 -2 3 -4 3 2 1];      % 输入信号
y = conv(h,x);              %h[n],x[n]都为有限长序列，且在线性时不变离散系统，求输出序列
n = 0:14;                   %输出序列的时间维度，9+7-1
stem(n,y);
xlabel('n'); ylabel('y[n]');
title('Output Obtained by Linear Convolution'); grid;

周期矩形脉冲信号的指数形式傅里叶级数

$$f(t)=\sum_{-\infty }^{\infty } F_{k} e^{jk\Omega _{0}t }$$ $$F_{k} =\int_{-\frac{T}{2} }^{\frac{T}{2} } f(t)e^{-jk\Omega _{0} t} dt=\frac{A\tau }{T} Sa(\frac{k\Omega _{0} \tau }{2} )$$ $$其中\Omega _{0} =\frac{2\pi }{T}$$

A=1; 
tau=0.1; %脉冲宽度
T=0.5;  %周期
Omega_0 = 2*pi/T;  %基波频率
K = 2*pi/tau/Omega_0; %%K=T/tau占空比的倒数，占空比：tau/T
k= 0:2*K;
F_k = A*tau/T*sinc(k*Omega_0*tau/2/pi);  % 系数F_k,Sa(x) = sin(x) / x,sinc(x) = Sa(pi * x)
F_k_mag=abs(F_k);      %F_k的幅度
F_k_phase = angle(F_k);  %F_k的相位
k= -2*K:2*K;
F_k_mag = [fliplr(F_k_mag(2:end)) F_k_mag];
F_k_phase = [ -fliplr(F_k_phase(2:end)) F_k_phase];
subplot(2,1,1)
stem(k*Omega_0, F_k_mag);
xlabel('k \Omega_o');
ylabel('magnitude');
grid
subplot(2,1,2)
stem(k*Omega_0, F_k_phase);
xlabel('k \Omega_o ');
ylabel('phase');
grid

前N次谐波叠加得到的周期矩形脉冲信号近似波形

吉布斯现象：将具有不连续点的周期函数（如矩形脉冲)进行傅立叶级数展开后，选取有限项进行合成。当选取的项数越多，在所合成的波形中出现的峰值越靠近原信号的不连续点。当选取的项数很大时，该峰值趋于一个常数

$$f_{N} (t)=\frac{A\tau }{T} +\sum_{k=1}^{N} \frac{2A\tau }{T}Sa(k\Omega _{0} \frac{\tau }{2} ) cosk\Omega _{0} t$$

t=-2:10^(-4):2; %采样间隔,会影响是否能采样到f_N的最大值,可增加采样点数可以减小误差
A=1;
tau=1; %脉冲宽度
T=2;  %周期
Omega_0 = 2*pi/T;  %基波频率
c0=A*tau/T; 
N=input('N=');
f_N=c0*ones(1,length(t));  %f_N(t)
for k=1:1:N  
     f_N = f_N + 2*A*tau/T*sinc(k*Omega_0*tau/2/pi) *cos(k*Omega_0*t);
end
plot(t, f_N);
xlabel('t');
ylabel('f_N(t)');
title(['N=',num2str(N)]) 
axis([-2 2 -0.2 1.2]);

对采样频率Fs的长度为N的序列x求频谱

function [X,f] = ctft1(x,Fs,N)
X=fftshift(fft(x,N))/Fs;
f=-Fs/2+(0:N-1)*Fs/N;
end

音频录制及存储

%------ 录制自己的语音
%创建一个保存音频信息的对象，包含采样率，时间和录制的音频信息。44100表示采样为44100Hz（可改为8000，11025，22050等，采样频率越高，录入的声音质量越好，相应需要的存储空间越大），16为用16bits存储，1为单声道，也可以改为2表示两通道）
R=audiorecorder(44100,16,1); %audiorecorder 创建一个 8000 Hz、8 位、1 声道的audiorecorder 对象。将返回对象的句柄。
record(R);%开始录制，在指令窗输入此指令后，对着麦克风说话，即可进行录制
disp('Start speaking.')
%录制者说一句话，说完后输入下面指令
pause(R);%暂停录制
play(R); % 播放录制的声音
%resume(R);%如果需要的话可以继续录制
%stop(R);%停止录制，录制结束

%------ 画出语音的时域波形
myspeech=getaudiodata(R);%得到刚录制的音频信号矢量,getaudiodata获取 audiorecorder 对象中录制的音频数据,将录制的音频数据作为双精度数组
返回
plot(myspeech);%画出语音波形
xlabel('时域样值'),ylabel('幅度'),title('语音波形');

%-----写入和读取声音文件
%myspeech.wav:音频文件完整的路径
audiowrite('myspeech.wav',myspeech,44100);%语音存储
%audiowrite（Y，Fs，NBITS，WAVEFILE） 将数据 Y写入由文件名 WAVEFILE 指定的windows wave文件。 windows wave文件具有采样率 FS Hz 和 NBITS 位数。 NBITS 必须为 8、16、24 或 32,Y=audioread（FILE） 读取由字符串 FILE 指定的 WAVE 文
件，以 Y 格式返回采样数据
[x,Fs]=audioread('myspeech.wav');

%-----计算并画出语音信号的幅度频谱
N=length(x);
[X,f]=ctft1(x,Fs,N);	
figure
plot(f,abs(X)/max(abs(X)));
xlabel('频率(Hz)'),ylabel('幅度谱'),title('语音幅度谱');
axis([0,5000,0,1]);grid    % 当只观察0~5000Hz内的频谱

DFT与IDFT

$$X[k]=\sum_{n=0}^{N-1} x[n]W_{N}^{kn} ,0\le k\le N-1 ,W_{N} =e^{-j\frac{2\pi }{N} }$$ $$x[n]=\frac{1}{N} \sum_{k=0}^{N-1} X[k]W_{N}^{-kn} ,0\le k\le N-1 ,W_{N} =e^{-j\frac{2\pi }{N} }$$

L=length(x);  
N=L;      % set samples number in frequency domain
w=2*pi/N*(0:N-1);    %  discrete frequency
%DFT
Xw=zeros(1,N);  % vector for storing DFT 
for k=1:N     
     Xw(k)=x*(exp(-j*w(k)*(0:L-1)'));  
end 

%IDFT
y=zeros(1,N);
for n=1:N
y(n)=(1/N)*Xw* (exp(j*(n-1)*w')); 
end

圆周移位

$$h[\left \langle n-m \right \rangle _{N} ]$$

function y=cirshift(x,M)   % M为序列圆周移位的样本数
%程序对一个有限长序列x圆周移位M个采样间隔得到新序列y
if abs(M)>length(x)
   M=rem(M,length(x));
end
if M<0
   M=M+length(x);
end
y=[x(M+1:length(x)) x(1:M)];

圆周卷积

$$y_{c} [n]=\sum_{m=0}^{N-1} g[m]h[\left \langle n-m \right \rangle _{N} ],0\le n\le N-1$$

function y=circonv(x1,x2);
L1=length(x1);
L2=length(x2);
if L1 ~=L2
  error('Sequences of unequal lengths')
end
y=zeros(1,L1);
x2tr=[x2(1) x2(L2:-1:2)];
for k=1:L1
   sh=cirshift(x2tr,1-k);
   h=x1.*sh;
   y(k)=sum(h);
end

圆周卷积求线性卷积（时域）

• 条件：圆周卷积的长度大于线性卷积长度

g1 = [1 2 3 4 5];g2 = [2 2 0 1 1];
g1e = [g1 zeros(1,length(g2)-1)];
g2e = [g2 zeros(1,length(g1)-1)];
ylin = circonv(g1e,g2e);
disp('Linear convolution via circular convolution = ');
disp(ylin);
y = conv(g1, g2);
disp('Direct linear convolution = ');
disp(y)

圆周卷积求线性卷积（频域，DFT）

g1 = [1 2 3 4 5];g2 = [2 2 0 1 1];
g1e = [g1 zeros(1,length(g2)-1)];
g2e = [g2 zeros(1,length(g1)-1)];
G1=fft(g1e);
G2=fft(g2e);
ylin=ifft(G1.*G2);
disp('Linear convolution via circular convolution = ');
disp(ylin);
y = conv(g1, g2);
disp('Direct linear convolution = ');
disp(y)

频率响应

系统的零状态响应的傅里叶变换与输入信号的傅里叶变换之比

H=freqs(b,a,w);
%b和a分别为分子多项式和分母多项式的系数向量；w为系统频率响应的频率范围的向量，如w1: dw: w2形式，w1为频率的起点，w2为频率的终点，dw为频率的采样间隔；H为系统频率响应的样值
[H,w]=freqs(b,a);
%输出从计算出的频率响应中自动选取200个频率点的频率响应的样值，其中b和a同上；H为保存200个频率点的系统频率响应的样值，w为保存200个频率点的位置
[H,w]=freqs(b,a,N);
输出从计算出的频率响应中自动选取N个频率点的频率响应的样值，其中b和a同上；N为频率点的个数；H为保存N个频率点的系统频率响应的样值，w为保存N个频率点的位置
freqs(b,a);
%以伯德图的方式绘制出系统的频率响应曲线，并不返回系统频率响应的样值

已知频率响应求冲激响应h(t)

$$h(t)=\frac{1}{2\pi } \int_{-\infty }^{\infty } H(j\Omega )e^{j\Omega t} d\Omega$$