【数字信号处理】

第一次

介绍Matlab的使用 略

第二次

Task1 求解差分方程

解差分方程$y_{n}-y_{n-1}+0.9y_{n-2}=x(n)$

$x(n)=\delta(n),x(n)=\mu(n),y(-2)=0$

求解方法——递推

差分方程：$\sum_{k=0}^{N} a_{k} y(n-k)=\sum_{k=0}^{M} b_{k} x(n-k)$

递归形式：$y(n)=\frac{1}{a_{0}}[\sum_{k=0}^{M} b_{k} x(n-k)-\sum_{k=1}^{N} a_{k} y(n-k)]$

对于本题：$y(n)=x(n)+y(n-1)-0.9y(n-2)$

​ 其中，$a=[1,-1,0.9],b=[1]$

验证方法——filter

编程实现

点击查看代码

clc;clear all
b=1;a=[1 -1 0.9];
N=100;
x=(-2:100);
h=zeros(1,N+2);
s=zeros(1,N+2);
%构造单位脉冲序列和单位阶跃序列
n=(-2:N);
delta=(n==0);
u=(n>=0);
%求单位脉冲响应
for k=3:N+3
    h(k)=b*delta(k)-a(2)*h(k-1)-a(3)*h(k-2);
end
subplot(2,1,1)
plot1 = plot(x,delta,'.');
hold on
plot2 = plot(x,h,'linewidth',1);
%y=impz(b,a,N);
%plot3 = stem((0:99),y','.');
y = filter(b,a,delta);
plot3 = stem(x,y,'.');
legend("单位冲激序列","单位冲激响应","通过impz验证")
title("求冲激响应")
hold off
%求单位阶跃响应
for k=3:N+3
    s(k)=b*u(k)-a(2)*s(k-1)-a(3)*s(k-2);
end
subplot(2,1,2);
plot4 = stem(x,u,'.');
hold on
plot5 = plot(x,s,'linewidth',1);
%y= stepz(b,a,N);
%plot6 = stem((0:99),y','.');
y = filter(b,a,u);
plot6 = stem(x,y,'.');
legend("单位阶跃序列","单位阶跃响应","通过stepz验证")
title("求阶跃响应")

Task2 自相关检测

雷达测距原理

编程实现

点击查看代码

clear all;clc;
%相位相关算法在延时检测中的应用
a=0.4;B=60;%衰减系数a,延时B
dt=0.2;
t=0:dt:100;
f=ft(6,t);%发射波
r=a*ft(6,t-B);%回波
L=length(t);
n=randn(1,L)/5;%噪声
rn=r+n;%带有噪声的回波
subplot(5,1,1);plot(t,f);ylabel('f(t)')
title('发出波')
subplot(5,1,2);plot(t,r);ylabel('r(t)')
title('回波')
subplot(5,1,3);plot(t,rn);ylabel('rn(t)')
title('叠加随机噪声的回波')
F=fftshift(fft(f));%求信号f的傅里叶变换F（jw）
R=fftshift(fft(r));%求信号x的傅里叶变换F（jw）,并将零频率点移至频谱的中心位置
G=conj(F).*R;%求F*（jw）R(jw)
Cp=ifft(G);%求傅里叶反变换
subplot(5,1,4);plot(t,abs(Cp));ylabel('Cp(t)')
title('相关算法在延时检测中的应用');
G=G./abs(G);%归一
Rn=fftshift(fft(rn));
G=conj(F).*Rn;G=G./abs(G);
Cp=ifft(G);
subplot(5,1,5);plot(t,abs(Cp));ylabel('Cpn(t)')
title('相位相关算法在延时检测中的应用');
function ft=ft(w,t)
    ft=heaviside(t)-heaviside(t-20);
    for k=1:w
        ft = ft+sin(pi*k/2*t).*(heaviside(t)-heaviside(t-20)) ;
    end
end

Task3 傅里叶变换

参考手动实现离散傅里叶正变换与逆变换(程序+例子) - 简书 (jianshu.com)

点击查看代码

data = [2 3 8 9 12 4 8 10];   % 原始离散时域信号
N = length(data);
 
% DFT需要的相关矩阵:
WN = exp(-i*2*pi/N);  % 常数
WN_nk = zeros(N)+WN;  % WN_kn矩阵(初始)
xk = data';     % 时域信号振幅(列矩阵)
E = zeros(N);   % 辅助的E(WN_kn的幂,单独拿出来算)
 
%%% 傅里叶正变换即结果 %%%
for row = 0:N-1
    for cow = 0:N-1
        E(row+1,cow+1) = row*cow;
        WN_nk(row+1,cow+1) = WN_nk(row+1,cow+1)^E(row+1,cow+1);
    end
end
 
Xk = WN_nk * xk;   % 频域结果        
figure;
subplot(1,2,1);stem(Xk,'b','fill');hold on
stem(fft(data),'r');    % Matlab自带语句,验证
legend('自己算的','Matalab自带语句验证')hold off
 
% 继DFT代码后跟着写
% IDFT需要的相关矩阵：
WN_nk_n = zeros(N)+WN;  % WN_kn矩阵(初始)
E_n = zeros(N);                   % 辅助的E(WN_kn的幂,单独拿出来算)
 
for row = 0:N-1
    for cow = 0:N-1
        E_n(row+1,cow+1) = -row*cow;     % 注意负号
        WN_nk_n(row+1,cow+1) = WN_nk_n(row+1,cow+1)^E_n(row+1,cow+1);
    end
end
 
xk_n = real((WN_nk_n * Xk)/N)';     % IDFT手写结果
subplot(1,2,2);stem(xk_n,'b','fill');hold on
stem(ifft(Xk),'r');                % Matlab自带语句,验证

第三次

Task1.计算DTFT

求$x(n)=(0.5)^{n}u(n)$的离散时间傅里叶变换

序列$x(n)$是绝对可加的，因此它的FT存在

$X(e^{jw})=\sum_{n=-\infty}^{\infty}x(n)e^{-jwn}=\sum_{n=0}^{\infty}(0.5)^{n}e^{-jwn}\\=\sum_{n=0}^{\infty}(0.5e^{-jw})^{n}=\frac{1}{1-0.5e^{-jw}}=\frac{e^{jw}}{e^{jw}-0.5}$
说明：

1.这是一个无限长的信号，但在Matlab中只能取N个有限制来近似

2.$w$是一个$(-\infty,\infty)$之间变化的实变量，在Matlab中尽量取间隔小一些模拟连续

3.周期性$X(e^{jw})=X(e^{j[w+2\pi]})$和对称性(对于实信号)$X(e^{jw})$共轭对称

其实只需要取$[0,\pi]$就可以知道$X$的全部信息，但为了观察我们选取区间$[\pi,\pi]$

点击查看代码

%要求：求x(n)=0.5^nu(n)的FT
%产生x(n)
clear;clc
N=500;
n=0:N;
xn=ones(1,N+1);
x=xn.*(power(0.5,n));
 
%画原始信号
figure;
subplot(3,2,[1,2]);plot(x);grid
title('Primary Singal');
K=1000;
k=-K:K;
 
%----------------FT变换-----------------
w=(2*pi/K)*k;
%1.方法一循环
X=ones(1,2*K+1);
for i=1:401
   X(i)=sum(x.*exp(-1i*w(i)*n)); 
end
%2.方法二矩阵
X_2=x*(exp(-1i*pi/N)).^(n'*k);
%3.公式手算验证
X_confirm=exp(1i*w)./(exp(1i*w)-0.5*ones(1,2*K+1));
T = (X==X_confirm);
T_2=(X_2==X_confirm);
 
%对结果分析展示
magX=abs(X);angX=angle(X);
realX=real(X);imagX=imag(X);
subplot(3,2,3);plot(k*2/K,magX);grid
title('|X(e^{-jw})|=|X(e^{jw})|(偶对称)');ylabel('Magnitude');
subplot(3,2,5);plot(k*2/K,angX);grid
title('∠X(e^{-jw}=-∠X(e^{-jw}(奇对称)');ylabel('Radians');
subplot(3,2,4);plot(k*2/K,realX);grid
title('Re[X(e^{-jw}]=Re[X(e^{jw}](偶对称)');ylabel('Real');
subplot(3,2,6);plot(k*2/K,imagX);grid
title('Im[X(e^{-jw}]=-Im[X(e^{-jw}](奇对称)');ylabel('Imaginary');
figure;
subplot(2,2,1);polarplot(w,magX);title('|X(e^{-jw})|=|X(e^{jw})|(偶对称)');
subplot(2,2,3);polarplot(w,angX);title('∠X(e^{-jw}=-∠X(e^{-jw}(奇对称)');
subplot(2,2,2);polarplot(w,realX);title('Re[X(e^{-jw}]=Re[X(e^{jw}](偶对称)');
subplot(2,2,4);polarplot(w,imagX);title('Im[X(e^{-jw}]=-Im[X(e^{-jw}](奇对称)');

Task2.计算DFT

计算$G_{4}(n)$的4/8/16点离散傅里叶变换

$X(K)=\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}kn}\quad (k=0,1,2,...,N-1)\\=\sum_{n=0}^{N-1}x(n)W_{N}^{nk},\quad W_{N}=e^{-j2\pi/N}$

以$N=4$为例

$X_{4}^{k}=\sum_{n=0}^{3}x(n)W_{k}^{nk};k=0,1,2,3;W_{4}=e^{-j2\pi/4}=-j$

$X_{4}^{0}=\sum_{n=0}^{3}(-j)^{0}=4$

$X_{4}^{1}=(-j)^{0}+(-j)^{1}+(-j)^{2}+(-j)^{3}=1+(-j)+(-1)+j=0$

$X_{4}^{2}=(-j)^{0}+(-j)^{2}+(-j)^{4}+(-j)^{-6}=1+(-1)+1+(-1)=0$

$X_{4}^{3}=(-j)^{0}+(-j)^{3}+(-j)^{6}+(-j)^{9}=1+j+-1-j=0$

说明：在Matlab里面因为精度的问题，0通常会是极小极小的数而不是真正的0，这就会导致用复数求相位不准确，所以在求相位之前需要先将复数按照某一误差归零，但是这一误差我还不知道怎么调好，所以16点4N点误归0了。

Task3.频谱-1

$x(n)=cos(0.48\pi n)+cos(0.52\pi n)$

1.取$0\leq n\leq 9$求$X_1(k)$频谱

2.将$x(n)$补90个0求$X_2(k)$频谱

3.取$0\leq n\leq 99$求$X_{3}(k)$

说明：

1.虽然Task2编的DFT虽然相频特性效果不佳，但这里只观察幅频特性，所以继续用$y*exp(-1i*n'*w)$，后来换回fft效果是一样的。

2.周期性$X(e^{jw})=X(e^{j[w+2\pi]})$和对称性(对于实信号)$X(e^{jw})$共轭对称。只需要取$[0,\pi]$就可以知道$X$的全部信息。

3.第一排图：样本太少，采样频率不足

3.第二排图：补零所得更长一些的DFT对原序列的离散时间傅里叶变换提供了更为密集的样本，从而给出了一种高密度的谱，但是它没有任何新的信息附加到这个信号，只是让第一排图更平滑。很明显$x(n)=cos(0.48\pi n)+cos(0.52\pi n)$的主要频率并不在$0.5\pi$而是在$0.48\pi$和$0.52\pi$。

4.第三排图：为了得到更高分辨率的谱，我们需要观察更多的数据，这样就可以较为准确地描述$x(n)=cos(0.48\pi n)+cos(0.52\pi n)$的幅频特性。

点击查看代码

%要求： x(n),0<=n<=9, 0<=n<=9+90,0<=n<=99 fft
clear;clc
n=0:1:99;x=cos(0.48*pi*n)+cos(0.52*pi*n);
%1.signal x(n),0<=n<=9 X_{1}(k)
n1=0:1:9;y1=x(1:1:10);
subplot(3,2,1);stem(n1,y1);title('x(n) 0<=n<=9');
k1=0:1:5;w1=2*pi/10*k1;
Y1=y1*exp(-1i*n1'*w1);magY1=abs(Y1(1:1:6));
subplot(3,2,2);stem(w1/pi,magY1);title('DTFT 幅频特性图[0,π]');
 
%2.Add zeros(1,90) X_{2}(k)
n2=0:1:99;y2=[x(1:1:10) zeros(1,90)];
subplot(3,2,3);stem(n2,y2);title('x(n) 0<=n<=9+90 zeros');
k2=0:1:50;w2=2*pi/100*k2;
Y2=y2*exp(-1i*n2'*w2);magY2=abs(Y2(1:1:51));
subplot(3,2,4);plot(w2/pi,magY2);title('DTFT 幅频特性图[0,π]');
 
%3.signal x(n),0<=n<=99 X_{3}(k)
subplot(3,2,5);stem(n,x);
title('x(n),0<=n<=99');xlabel('n')
k=0:1:50;w=2*pi/100*k;
X=x*exp(-1i*n'*w);magX=abs(X(1:1:51));
subplot(3,2,6);plot(w/pi,magX);title('DTFT 幅频特性图[0,π]');

Task4.频谱-2

模拟信号$x_{a}=cos(200\pi t)+sin(100\pi t)$采样率$f_{s}=400Hz$采0.04s/0.16s/0.32s求频谱

采样周期$T_{s}=1/f_{s}$

以相距$T_{s}$秒的采样间隔对$x_{a}$采样，得到离散时间信号$x(n)=x_{a}(nT_{s})$
随着采样点数的增多，得到更高分辨率的谱，这样就可以更为准确地描述信号的幅频特性。

点击查看代码

%要求：fs=400,t=0.04,0.16,0.32 dft
clear;clc
fs=400;
Ts=1/fs;
t1=0.04;N1=t1*fs;n1=0:1:N1;x1=cos(200*pi*n1*Ts)+sin(100*pi*n1*Ts);
t2=0.16;N2=t2*fs;n2=0:1:N2;x2=cos(200*pi*n2*Ts)+sin(100*pi*n2*Ts);
t3=0.32;N3=t3*fs;n3=0:1:N3;x3=cos(200*pi*n3*Ts)+sin(100*pi*n3*Ts);
%DTFT
K=500;k=0:1:K;w=pi*k/K;
X1=x1*exp(-1i*n1'*w);X1=real(X1);
X2=x2*exp(-1i*n2'*w);X2=real(X2);
X3=x3*exp(-1i*n3'*w);X3=real(X3);
%展示
subplot(3,2,1);stem(n1,x1);title('f=400Hz,t_{1}=0.04')
subplot(3,2,2);plot(w/pi,X1);title('X_{1}(k)')
subplot(3,2,3);stem(n2,x2);title('t_{2}=0.16')
subplot(3,2,4);plot(w/pi,X2);title('X_{2}(k)')
subplot(3,2,5);stem(n3,x3);title('t_{3}=0.32')
subplot(3,2,6);plot(w/pi,X3);title('X_{3}(k)')

第四次

模仿钢琴音色

点击查看代码

clear;clc
%made by tty 2021.11.15
%泛频参考：wuli玩 | MATLAB的音乐物理学（上）https://mp.weixin.qq.com/s/9XQ1mzeGzE0-jL9ijoIcww
%《钢琴及其调律》廖志坚 附录
%% Musical Offering Canon A2 BWV1079
music.pitch=[
    "C" "#D"...
    "G" "#G"...
    "B1" "G"...
    "#F" "F"...
    "E" "#D"...
    "D" "#C" "C" "B1"...
    "G1" "C" "F" "E"... 
    "#D" "D"...
    "C" "#D"...
    "G" "F" "G" "c" "G" "#D" "D" "#D"...
    "F" "G" "A" "B" "c" "#D" "F" "G"...
    "#G" "D" "#D" "F" "G" "F" "#D" "D"...
    "#D" "F" "G" "#G" "#A" "#G" "G" "F"...
    "G" "#G" "#A" "c" "#c" "#A" "#G" "G"...
    "A" "B" "c" "d" "#d" "c" "B" "A"...
    "B" "c" "d" "#d" "f" "d" "G" "d"...
    "c" "d" "#d" "f" "#d" "d" "c" "B"...
    "c" "G" "#D" "C" "Sh" ...
    ];
 
music.Rhythm=[
    2 2 ...
    2 2 ...
    4/3 4 ...
    4/3 4 ...
    4/3 4 ...
    4 4 4 4 ...
    4 4 4 4 ...
    2 2 ...
    2 2 ...
    8 8 8 8 8 8 8 8 ...
    8 8 8 8 8 8 8 8 ...   
    8 8 8 8 8 8 8 8 ... 
    8 8 8 8 8 8 8 8 ...   
    8 8 8 8 8 8 8 8 ...   
    8 8 8 8 8 8 8 8 ... 
    8 8 8 8 8 8 8 8 ...   
    8 8 8 8 8 8 8 8 ... 
    8 8 8 8 2 ...
    ];
 
%% 播放
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%设置参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs=44100;%mp3标准采样品频率表44100
%四分音符为一拍，一拍0.5s,这个曲子4/4拍，拍的时长为2s
Duration=2;
T=0:Fs^-1:Duration;
L=length(T);
%%%%%%%%%%%%%%%%%%%%%%%%%%%分两个声部（左右手）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
right=zeros(1, L);left=zeros(1, L);
right=piece(Duration,music,Fs,2);                   %右手高音，频率高一些2*f
left=piece(Duration,music,Fs,1);                      %右手低音，频率低一些f
right = smoothdata(right); %做简单的平滑处理，希望没有剧烈变化又不希望旋律失真 
left = smoothdata(left); 
pplay = fliplr(left);                 %根据螃蟹卡农的特点，正放和倒放都成旋律
stereo(:,2)=right;                                   %将左右手旋律存入双声道
stereo(1:length(right),1)=pplay;
%ppplay=play+pplay;%error
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%写入文件%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sound(stereo,Fs);
% audiowrite("D:\大三上\3.5数字信号处理\钢琴\正向.wav",right,Fs);
% audiowrite("D:\大三上\3.5数字信号处理\钢琴\反向.wav",left,Fs);
% audiowrite("D:\大三上\3.5数字信号处理\钢琴\螃蟹卡农.wav",stereo,Fs);
 
%根据键盘名获取对应的频率
function frequency=hz(name)
%音列
KeyBoard.name=[
    "A2" "#A2" "B2"...                                                     %大字二组
    "C1" "#C1" "D1" "#D1" "E1" "F1" "#F1" "G1" "#G1" "A1" "#A1" "B1"...    %大字一组
    "C" "#C" "D" "#D" "E" "F" "#F" "G" "#G" "A" "#A" "B"...                %大字组(1-27为低音区)
    "c" "#c" "d" "#d" "e" "f" "#f" "g" "#g" "a" "#a" "b"...                %小字组
    "c1" "#c1" "d1" "#d1" "e1" "f1" "#f1" "g1" "#g1" "a1" "#a1" "b1"...    %小字一组
    "c2" "#c2" "d2" "#d2" "e2" "f2" "#f2" "g2" "#g2" "a2" "#a2" "b2"...    %小字二组（28-63为中音区）
    "c3" "#c3" "d3" "#d3" "e3" "f3" "#f3" "g3" "#g3" "a3" "#a3" "b3"...    %小字三组
    "c4" "#c4" "d4" "#d4" "e4" "f4" "#f4" "g4" "#g4" "a4" "#a4" "b4"...    %小字四组
    "c5"...                                                                %小字五组（64-88为高音区）
    ];
KeyBoard.rest=("Sh");
%调音
if (ismember(KeyBoard.rest,name))
    frequency=0;
else
    n=find(KeyBoard.name==name);                                           %根据谱号找其在整个键盘中的位置
    c=find(KeyBoard.name=="c");%28
    %a=find(KeyBoard.name=="a");sigma=2^(1/12);   
    pure=[261.63 275.59 290.63 310.05 326.97 348.83 367.86 391.99 413.39 435.95 465.09 490.47];   %纯律系数
    tw=[261.63 277.18 293.66 311.13 329.63 349.23 369.99 391.99 415.30 440.00 456.16 493.88];     %用十二平均律系数                                                                           
    fifths=[261.63 279.11 294.34 314.02 331.15 349.23 372.37 392.45 418.68 441.55 471.04 496.74]; %五度相生律系数,
    if n < 28
        p=mod((n-c),12)+1;
        group=(n-p-2)/12-1;
        frequency=pure(p)*2^(group);
    elseif n < 64
        %frequency=440*sigma^(n-a);
        p=mod((n-c),12)+1;
        group=(n-p-2)/12-1;
        frequency=tw(p)*2^(group);
    else 
        p=mod((n-c),12)+1;
        group=(n-p+2)/12-1;
        frequency=fifths(p)*2^(group);
    end
end
end
 
%合成音符
function part=note(Fs,f,T)
    t=0:Fs^-1:T;
    L=length(t);
    %不同频率范围的泛音列
    if Fs<=33
        lvl_wave = [10000 2190 6095 2297 2131 1139 1745 576.8 724.8 883.4,...
            207.9 355.4 403.9 146.2 901 392.6 212.5 192.4 368.8 234.9 99.55 198.8 393.8 72.65 89.79] / 10000;
    elseif Fs>33&&Fs<=524
        lvl_wave = [10000 3007 6964 587.5 1185 660.4 1206 370.9 580.3 501.2 ...
            260 224.3 57.98 319.9 685.9 550.4 200.3 200.3 244 221.4 115.9 71.42 299.3] / 10000;
    elseif Fs>524&&Fs<=1047
        lvl_wave = [10000 6447 2793 3265 938.3 1371 1284 1.909 2549 1075 ...
            44.95 1211 690.5 1085 566.8 428.9] / 10000;
    elseif Fs>1047&&Fs<=2094
        lvl_wave = [10000, 2371, 885.3, 258.6] / 100000;
    else
        lvl_wave = [10000, 541] / 100000;
    end
    %依据泛音列合成泛音
    overtone = zeros(1, L);                                                %泛音初始化
    for i = 1 : length(lvl_wave)
        overtone = overtone  + lvl_wave(i) * sin(pi*f*i*t);                %模仿泛音
    end
    overtone=6*overtone;
    %包络
    envelope1 =(1 - exp(-80*(t))) ./ (1 + exp(-80*(t)));                   % 包络函数1，指数渐入的效果
    envelope1= [(zeros(fix(L/100)-1,1))' envelope1((fix(L/100)):end)];     % 模仿wav前面有一些空白
    envelope2 = (t.^0.01.*exp(-2*t));                                      % 包络函数2，指数渐出的效果
    part=overtone.*envelope1.*envelope2; 
end
 
%合成一个声部
function single=piece(Duration,music,Fs,m)
single=[];
for i=1:length(music.pitch)
    f=hz(music.pitch(i));%基音频率，单位赫兹
    t=Duration*(1/music.Rhythm(i));%每个音对应的音长
    part=note(Fs,m*f,t);
    single=[single,part];
end
end
 
%% 其它编程
% %根据键盘名获取对应的频率
% function frequency=hz(name)
% %音列
% KeyBoard.name=[
%     "A2" "#A2" "B2"...                                                     %大字二组
%     "C1" "#C1" "D1" "#D1" "E1" "F1" "#F1" "G1" "#G1" "A1" "#A1" "B1"...    %大字一组
%     "C" "#C" "D" "#D" "E" "F" "#F" "G" "#G" "A" "#A" "B"...                %大字组(1-27为低音区)
%     "c" "#c" "d" "#d" "e" "f" "#f" "g" "#g" "a" "#a" "b"...                %小字组
%     "c1" "#c1" "d1" "#d1" "e1" "f1" "#f1" "g1" "#g1" "a1" "#a1" "b1"...    %小字一组
%     "c2" "#c2" "d2" "#d2" "e2" "f2" "#f2" "g2" "#g2" "a2" "#a2" "b2"...    %小字二组（28-63为中音区）
%     "c3" "#c3" "d3" "#d3" "e3" "f3" "#f3" "g3" "#g3" "a3" "#a3" "b3"...    %小字三组
%     "c4" "#c4" "d4" "#d4" "e4" "f4" "#f4" "g4" "#g4" "a4" "#a4" "b4"...    %小字四组
%     "c5"...                                                                %小字五组（64-88为高音区）
%     ];
% KeyBoard.rest=("Sh");
% %调音
% if (ismember(KeyBoard.rest,name))
%     frequency=0;
% else
%     n=find(KeyBoard.name==name);                                               %根据谱号找其在整个键盘中的位置
%     a=find(KeyBoard.name=="a");sigma=2^(1/12); 
%     frequency=440*sigma^(n-a);
% end
% end
% music.pitch=["C" "C" "G" "G" "A" "A" "G" "Sh" "F" "F" "E" "E" "D" "D" "C"];
% music.Rhythm=[4 4 4 4 4 4 4 8 4 4 4 4 4 4 4];
% music.pitch=[
%     "E" "G" "G" "A" "G" "E" "C" "D" "E" "E" "D" "C" "D" ...
%     "E" "G" "G" "A" "G" "E" "C" "D" "E" "E" "D" "D" "C" ...
%     ];
% music.Rhythm=[
%     4 4 2 8 ...
%     4 4 2 8 ...
%     4 4 4 4 ...
%     1 ...
%     4 4 2 8 ...
%     4 4 2 8 ...
%     4 4 4 4 ...
%     1 ...
%     ];

音乐：螃蟹卡农（点击下载）

第五次

实序列的FFT算法

合二为一：用一个N点FFT计算两个实序列的FFT

$\begin{aligned} x_{1}(n)+jx_{2}(n)\leftrightarrow X(k)\\X_{1}(k)=X_{ep}=\frac{1}{2}[X(k)+X^{*}(-k)]=\frac{1}{2}[X(k)+X^{*}(N-k)]\\X_{2}(k)=-jX_{op}=-j\frac{1}{2}[X(k)-X^{*}(-k)]=-j\frac{1}{2}[X(k)-X^{*}(N-k)]\\(k=0,1,...,N-1)\end{aligned}$

一个小例子：

设两个长度$N=4$的实序列：

$\begin{aligned}x_{1}(n)_{n=0}^{3}=(1,1,3,1)\\x_{2}(n)_{n=0}^{3}=(2,1,2,2)\end{aligned}$

解：令$x(n)=x_{1}(n)+x_{2}(n)(n=0,1,2,3)$计算格式如下图

点击查看代码

clear;clc;close all;    %初始化
%实序列的FFT算法对比
 
%参数初始化设置
N = 4;                    %FFT点数
%fft_time1=zeros(1,N);     %存储运行所需时间的向量，1指“合二为一”方法
%fft_time2=zeros(1,N);     %存储运行所需时间的向量，2指“直接FFT”方法
x1 = [1,1,3,1];           %第1个实序列
x2 = [2,1,2,2];           %第2个实序列
 
%合二为一：用一个N点FFT计算两个实序列的FFT
tic
x = complex(x1, x2);      %根据x1、x2构造复数序列
X=fft(x);                 %对复序列做FFT
Xf=flip(X);               %为了求X*(N-k)，先做翻转
Xf=[Xf(end) Xf(1:end-1)]; %X(0)即X(N),坑
Xc=conj(Xf);              %取复共轭
X1=0.5*(X+Xc);            %按照公式，时域的实部对应频域的对称部分
X2=-1j*0.5.*(X-Xc);       %时域的虚部对应频域的反对称部分
toc
 
%与分别直接做FFT对比验证
tic
X11=fft(x1);              
X22=fft(x2);
toc
confire_1=(X1==X11);
confire_2=(X2==X22);

点击查看代码

clear;clc;close all;      %初始化
%实序列的FFT算法对比
 
%参数初始化设置
Nmax = 2048*2;               %FFT点数从1变到2048*2看两种方法用时变化
fft_time1=zeros(1,Nmax);     %存储运行所需时间的向量，1指“直接FFT”方法
fft_time2=zeros(1,Nmax);     %存储运行所需时间的向量，2指“合二为一”方法
 
for n=1:1:Nmax
    x1 = rand(1,n);
    x2 = rand(1,n);
    t1=tic;Realfft(x1,x2);fft_time1(n)=toc(t1);
    t2=tic;Myfft(x1,x2);fft_time2(n)=toc(t2);
end
n=1:1:Nmax;
plot(n,fft_time1,':b');
hold on;
plot(n,fft_time2,':r');
legend('Realfft','Myfft')
xlabel('N');ylabel('Time in Sec.')
title('FFT execution times')
 
%“直接FFT”方法
function [X1,X2]=Realfft(x1,x2)
    X1=fft(x1);              
    X2=fft(x2);
end
%“合二为一”方法
function [X1,X2]=Myfft(x1,x2)
    x = complex(x1, x2);      %根据x1、x2构造复数序列
    X=fft(x);                 %对复序列做FFT
    Xf=flip(X);               %为了求X*(N-k)，先做翻转
    Xf=[Xf(end) Xf(1:end-1)]; %X(0)即X(N),坑
    Xc=conj(Xf);              %取复共轭
    X1=0.5*(X+Xc);            %按照公式，时域的实部对应频域的对称部分
    X2=-1j*0.5.*(X-Xc);       %时域的虚部对应频域的反对称部分
end

第六次

Butterworth Filter

一个N阶低通滤波器的幅度平方响应：

$$|H_{c}(j\Omega)|^{2}=\frac{1}{1+(\frac{j\Omega}{j\Omega_{c}})^{2N}}$$

$\Omega_{c}$表示截止频率（3dB截止频率）

阶数N的大小主要影响通带幅频特性的平坦程度和过渡带、阻带的下降速度，它由技术指标$\Omega_{p},\Omega_{s},\alpha_{p}$和$\alpha_{s}$

技术指标
$\Omega_{p}$	通带截止频率
$\Omega_{s}$	阻带截止频率
$\alpha_{p}$	通带最大衰减（即通带$[0,\Omega_{p}]$中允许$A(\Omega)$的最大值）	dB
$\alpha_{s}$	阻带最小衰减（即阻带$\Omega\geq \Omega_{s}$上允许$A(\Omega)$的最小值）	dB

点击查看代码

clear;clc;
 
%对输入的两个正弦信号x1(t)、x2(t)进行采样
f1=200;     %正弦信号x1(t)的频率
f2=1700;    %正弦信号x2(t)的频率
Fs=10000;           %采样频率Fs
Ts=1/Fs;            %采样间隔Ts
t=0:Ts:2;           %“录音”2s
x1=cos(2*pi*f1*t);  %数字信号x1
x2=cos(2*pi*f2*t);  %数字信号x2
x=x1+x2;
L=length(t);
Y1 = fft(x);
P2 = abs(Y1/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
subplot(4,1,1),plot(f/1000,P1) 
title('Single-Sided Amplitude Spectrum of x(t)');xlabel('f/kHz');ylabel('|P1(f)|');xlim([0,2])
 
%设置低通滤波器指标
%通带截止频率 %阻带截止频率 %通带最大衰减 %阻带最大衰减
wp=2*pi*600;ws=2*pi*800;Rp=3;As=30;
[N,wc]=buttord(wp,ws,Rp,As,'s');
[z,p,k]=buttap(N);
[bp,ap]=zp2tf(z,p,k);
[B,A]=lp2lp(bp,ap,wp);
% [B,A]=butter(N,wc,'s');
k=0:L/2;wk=2*pi*f;
Hk=freqs(B,A,wk);
subplot(4,1,2),plot(f/1000,20*log10(abs(Hk)));grid on;hold on
pha=angle(Hk)*180/pi;plot(f/1000,20*log10(abs(pha)));
title(sprintf('N = %d Butterworth Lowpass Filter',N));xlabel('f/kHz');ylabel('-A(f)/dB');xlim([0,2])
 
%除去高频信号
A=2*abs(Y1(1:L/2+1)).*(abs(Hk));
subplot(4,1,3),plot(f/1000,A);xlim([0,2]);grid on
%还原信号
A=[A(1:end-1),flip(A)];
z2=real(ifft(A));
subplot(4,2,7),plot(f(1:200)/1000,z2(1:200)*L/2)
subplot(4,2,8),plot(f(1:200)/1000,x1(1:200)*L/2)

第七次

FIR实现线性相位

FIR实现非线性相位

点击查看代码

clear all;clc;
%设计滤波器
%Design a 31-tap linear-phase lowpass filter. Display its magnitude and phase responses.
b = cfirpm(30,[-1 -0.5 -0.4 0.7 0.8 1],@lowpass);
fvtool(b,1,'OverlayedAnalysis','phase')
%Design a nonlinear-phase allpass FIR filter of order 22 with frequency response given approximately by an efunction
% n = 22;                              % Filter order
% f = [-1 1];                          % Frequency band edges
% w = [1 1];                           % Weights for optimization
% b = cfirpm(n,f,'allpass',w,'real');  % Approximation
 
[audio,fs]=audioread('D:\大三上\3.5数字信号处理\上机matlab\第七次\boliping.0033.wav');%声音读取
audio = audio(:,1); %双通道变单通道
T = 1/fs;%采样间隔
n = length(audio);
t = (0:n-1)*T;%时间轴
f = (0:n-1)/n*fs;%频率轴
k=-fs+1:fs;
w=(2*pi/fs)*k;
%快速傅里叶变换
audio_fft = fftshift(fft(audio,n)); 
audioab = abs(audio_fft);
subplot(2,2,1)
plot(k/fs,20*log10(audioab))  %幅值变换为分贝单位
title('fft后信号幅值');ylabel 'A / dB';ylim([-100,50])
subplot(2,2,2)
plot(k/fs,angle(audio_fft)/pi)
title('fft后信号相位');ylabel 'Phase / \pi'
%滤波
z=filter(b,1,audio);
z_fft= fftshift(fft(z));     %滤波后的信号频谱
audiowrite('6.wav',abs(z),fs)
hold on;
subplot(2,2,3)
plot(k/fs,20*log10(z_fft))  %幅值变换为分贝单位
title('滤波后fft信号幅值');ylabel 'A / dB';ylim([-100,50])
subplot(2,2,4)
plot(k/fs,angle(z_fft)/pi)
title('滤波后fft信号相位');ylabel 'Phase / \pi'