【数字信号处理】
第一次
介绍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'
