omega=logspace(-3,1,500);  %logspace函数生成对数间隔向量
h=1./(1+1i*2*omega);
figure;
subplot(2,1,1);grid on;
semilogx(omega,20*log10(abs(h)));%在x轴取对数
title('对数幅频特性');
subplot(2,1,2);grid on;
semilogx(omega,angle(h));
title('相频特性');

2.试画出频响函数

的对数幅频特性。

w=0:0.02:100;
magHw=abs(3./(3*1i*w+2-w.^2));
semilogx(w,magHw);
title('对数幅频特性')
xlabel('Frequency in rad/sec-log scale');
ylabel('Magnitude of Vout/Vin');
grid

3.已知信号为

，用MATLAB编程实现该信号经冲激脉冲，抽样得到的抽样信号fs(t)及其频谱。令参数E=5，τ=0.5，采用抽样间隔

s=1;
dt=0.1;
Ts=0.2;
t1=-4:dt:4;
ft=(5*(1+cos(2*pi*t1))/2).*(heaviside(t1+0.5)-heaviside(t1-0.5));%阶跃函数的表达
subplot(221);plot(t1,ft),
grid on;
axis([-2 2 -0.1 5.1]);
xlabel('Time(sec)'),ylabel('f(t)');
title('信号ft');
N=500;
k=-N:N;
W=pi*k/(N*dt);
Fw=dt*ft*exp(-1i*t1'*W);
subplot(222);plot(W,abs(Fw)),
grid on;axis([-30 30 -0.1 1.1*pi]);
xlabel('\omerga'),ylabel('F(w)');
title('信号ft的频谱');
Ts=0.2;
t2=-4:Ts:4;
fst=(5*(1+cos(2*pi*t2))/2).*(heaviside(t2+0.5)-heaviside(t2-0.5));
subplot(223);plot(t1,ft,':'),hold on;
stem(t2,fst),grid on;
axis([-2 2 -0.1 5.1]);
xlabel('Time(sec)'),ylabel('fs(t)');
title('抽样后的信号'),hold off;
Fsw=Ts*fst*exp(-1i*t2'*W);
subplot(224);plot(W,abs(Fsw)),grid on;
axis([-30 30 -0.2 1.1*pi]);
xlabel('\omega'),ylabel('Fs(w)');
title('抽样后的信号的频谱')

4.对题3获得的抽样信号，采用截止频率为4pi的低通滤波器对其滤波后重建信号f(t)，并计算重建信号与原升余弦脉冲信号的绝对误差。

E=5;
wc=4*pi;
Ts=0.2;
n=-100:100;
nTs=n*Ts;
fs=(E*(1+cos(2*pi*nTs))/2).*(heaviside(nTs+0.5)-heaviside(nTs-0.5));
t=-4:0.1:4;
ft=fs*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));
t1=-4:0.1:4;
f1=(5*(1+cos(2*pi*t1))/2).*(heaviside(t1+0.5)-heaviside(t1-0.5));
subplot(311);
plot(t1,ft,':'),hold on;
stem(nTs,fs),grid on;
axis([-2 2 -0.1 5]);
xlabel('nTs'),ylabel('f(nTs)');
title('抽样间隔Ts时的抽样信号f(nTs)');
hold off;
subplot(312);
plot(t,ft),grid on;
axis([-4 4 -0.1 5]);
xlabel('t'),ylabel('f(t)');
title('由f(nTs)信号重建得到原信号');
error=abs(ft-f1);subplot(313);
plot(t,error),grid on;
xlabel('t'),ylabel('error(t)');
title('重建信号与原信号的绝对误差')

5.已知调制信号f(t)=cos5πt，载波信号fc(t)=cos60πt，编程画出调制与解调过程中的波形图与频谱图。

clear clc
fs=40;
Fs=400;
Fc=30;
N=400;
n=0:N;
t=n/Fs;
xt=cos(pi*5*t);
xct=cos(2*pi*Fc*t);
yt=xt.*xct;
Xw=fftshift(abs(fft(xt,512)));
Yw=fftshift(abs(fft(yt,512)));
ww=-256:255;
ww=ww*Fs/512;
figure,subplot(2,1,1),plot(t,xt),title('调制信号波形');
subplot(2,1,2),plot(ww,Xw),title('调制信号频谱');
figure,subplot(2,1,1),plot(t,yt),title('已调信号波形');
subplot(2,1,2),plot(ww,Yw),title('已调信号频谱');
ylt=yt.*xct;
figure,subplot(2,1,1),plot(t,ylt),title('解调过程中间信号波形');
Ylw=fftshift(abs(fft(ylt,512)));
subplot(2,1,2),plot(ww,Ylw),title('解调过程中间信号频谱');
h=fir1(20,40/200,hamming(21));
figure,freqz(h,1),title('滤波器频率特征');
y2t=filter(h,1,ylt);
Y2w=fftshift(abs(fft(y2t,512)));
figure,subplot(2,1,1),plot(t,y2t),title('解调信号的波形');
subplot(2,1,2),plot(ww,Y2w),title('解调信号的频谱')

%运行程序，观察结果，认真完成以下问题。

% ①如何理解任何周期信号都可以分解成正弦信号加权和的形式；

%答：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。傅里叶级数就是加权的权重，傅里叶级数分解的方法可以把任意周期信号表示为正弦和余弦信号的和，余弦可以表示为相移了90°的正弦。通过实验也可以很明显地看出来，最大谐波次数为400时，信号趋近于周期方波信号。

% ②怎样理解非周期信号的傅立叶变换可以由周期信号的傅立叶系数被其基频归一化获得；

%答：非周期信号的傅里叶变换可以看作周期为无穷大的周期信号，此时频谱间隔2*pi/w趋于零，谱线变为连续的。为了比较无穷小之间的比例，引入归一化的概念，将傅里叶系数基频归一化在图像处中，归一化后可以对频域能量图像更好的量化(将频域能量重新量化到0-255范围的灰度图像)

% ③合理解释信号的高频信息对应其细节，信号的低频信息对应其概貌；

%答：二维的图像可以分解成不同的频率成分。其中，低频成分描述大范围的信息，而高频成分描述具体的细节。在灰度图像中，亮度变化小的区域主要是低频成分，而亮度变化剧烈的区域 (比如物体的边缘)主要是高频成分。而低频信号对应信息中变化慢的那一部分，高频信号对应信息中变化慢的那一部分，概貌轮廓的变化比较慢，细节处的变化比较快，因此信号的高频信息对应其细节，信号的低频信息对应其概貌。

% ④解释时域平移信号对其频域的幅度谱不产生影响的原因；

%答：信号在时域左移/右移 t0,则频域信号会乘exp(±jwt0)，它的幅值为1，即不影响其幅度谱。

% ⑤叙述音频信号在时域进行尺度变换后，对听觉产生的影响，并解释其原因；

%答：音频信号在时域进行尺度扩展后，声音变得缓和，尺度压缩后，声音变得短促。这会影响音频的响度，具体表现为扩展后音量变小，压缩后音量变大

% ⑥根据音频信号在时域微分信号的频谱图，解释听到的该声音与原声音（微分之前）不同的原因；

%答：对信号进行微分后，频域信号相当于乘以jw，高频信号幅度增加，低频信号幅度减小。所以微分后信号中的高频分量更多，低频分量更少，起到一定的锐化作用。

% ⑦声音信号在频域平移1Hz后，根据其傅立叶逆变换模信号对应的时域图，解释听到声音的变化，并在理论上解释其原因；

%答：声音的大小会有一定的波动，因为频域平移，相当于时域信号乘以exp(jwt),在实数域上表现为乘以三角函数，其大小会发生周期性波动。

% ⑧讨论相位失真对人的听觉和视觉产生的影响；

%答：视觉上表现为看到的二维图像的边缘锐化减小，图像更加模糊，不同灰度值之间的区别也有所减小，听觉上表现为不同频率声音的区别变小

% ⑨声音信号重构成二维图像信号后，讨论乐音和噪声对视觉效果的不同；

%答：乐声信号中的高频信号多，对应的图像细节较多，噪声的低频信号更多，图像更加平滑。

% ⑩图像信号展开成一维信号后，讨论有规律的纹理图像、任意图像和噪声图像对应的一维音频信号对听觉效果的不同。

%答：有规律的纹理图像转化为的音频信号中规律部分较多，或许可以组成某种固定的和旋等重复部分，让人印象深刻，而任意图像这种随机性就增加，规律感就会减少，而噪声图像由于有较多的杂乱部分，听觉效果较差。

% 7.1.2 分别采用前 4、40、400 项，画出周期矩形脉冲信号的近似图

clear;close all; format compact
t=-1:0.001:1;
E=1;T=1;w=2*pi/T;
ft=square(2*pi*t,50);           %写出函数表达式
figure(1);
h=figure(1);
set(h,'name','7.1.2 分别采用前 4、40、400 项，画出周期矩形脉冲信号的近似图');
set(h,'position', [0.2 0.3 0.6 0.6],'Units','normalized' );%以相对大小的方式绘图
set(h,'position', [0.2 0.3 0.6 0.6],'Units','normalized' );
subplot(2,2,1),plot(t,ft),grid on
axis([-1,1,-1.5,1.5])%设置坐标轴 
title('周期方波信号')
n_max=[4 40 400];
N=length(n_max);
for k=1:N   %通过for循环的方式绘图
    n=1:2:n_max(k);
    b=4*E./(pi*n);  %设置幅度
    x=b*sin(2*pi*n'*t); %设置函数表达式
    subplot(2,2,k+1),plot(t,x),grid on
    axis([-1,1,-1.5,1.5])
    title(['最大谐波次数=',num2str(n_max(k))])
end  %结束循环标志
pause;

%7.1.4 任意周期信号的周期核函数傅立叶系数的计算及重构算法时间渐进合成演示

Um=0.1;T=1;w=2*pi/T;  %分别对最大值，周期和频率进行初始化
num_noise=100;   %定义噪声
%N=input('取的谐波次数 N= ');
N=60;   %设定前N项
Fs=1000;
t=0:1/Fs:T;
num_points=numel(t); %设置时间点的个数
dt=T/(num_points-1);  %设置时间点间隔
un=Um*(sin(1*2*pi*t));  %分解的波形
n=randn(num_noise-2,1);  %随机生成矩阵
n=[0;n;0];
n = interp1([0:num_noise-1],n,linspace(0,num_noise-2,num_points),'linear'); %此处采用线性插值法插入前面的噪声信号
un=un+0.8*n;
figure(2);
h=figure(2);
set(h,'position', [0.2 0.05 0.6 0.8],'Units','normalized' ); %以相对尺寸的方式画图
set(h,'position', [0.2 0.05 0.6 0.8],'Units','normalized' );
subplot(3,1,1);plot(t,un,'LineWidth',3);
axis([-0.05*T  T+0.05*T  -1.8  1.8]);
grid on;
hold on;

a0=trapz(un)*dt/T*1;  %初始化参数a0、b0、Phi0
b0=0; 
A0=abs(a0);
Phi0=0;
for k=1:N   %以循环的方式求解a(k) 、b(k)、 A(k)、 Phi(k)
    a(k)=trapz(un.*cos(k*w*t))*dt/T*2;
    b(k)=trapz(un.*sin(k*w*t))*dt/T*2;
    A(k)=sqrt(a(k)^2+b(k)^2);
    Phi(k)=atan2(-b(k),(a(k)+eps));
end
subplot(3,1,2);stem(0,A0);  axis([0 N -inf,inf]);title('单边幅度谱')
subplot(3,1,3);stem(0,Phi0); axis([0 N -inf,inf]);title('单边相位谱')
Reu=A0;
subplot(3,1,1);plot(t,Reu+0*t,'g');
for k=1:N
    clf;
    subplot(3,1,1);plot(t,un,'LineWidth',3);grid on;hold on;
    axis([-0.01*T  T+0.01*T  -1.8  1.8]);  %初始化坐标轴
    title(strcat('原图和前',num2str(k),'次谐波之和'));
    subplot(3,1,2);stem(0:k,[A0,A(1:k)]);title('单边幅度谱');axis([0 N -inf,inf]);grid on;
    subplot(3,1,3);stem(0:k,[Phi0,Phi(1:k)]);title('单边相位谱');axis([0 N -3.2 3.2]);grid on;
    ad=A(k)*cos(k*w*t+Phi(k));
    subplot(3,1,1);plot(t,ad,'k');
    subplot(3,1,1);plot(t,Reu,'r','LineWidth',2);
    if k==18
        pause;
    end
    pause(0.1);
    Reu=Reu+ad;  
end
Psll=trapz(un.^2)*dt/T;     % 原信号能量
Ps12=A0^2+sum(A(1:end).^2/2);  % 前N项谐波能量
pause;

%7.2.2 比较非周期信号傅立叶变换与周期信号傅立叶系数的关系

figure(3);
h=figure(3);
set(h,'name','7.1.4 任意周期信号的周期核函数傅立叶系数的计算及重构算法时间渐进合成演示');
set(h,'position', [0.2 0.05 0.6 0.4],'Units','normalized' );  %以相对尺寸的方式绘图
set(h,'position', [0.2 0.05 0.6 0.4],'Units','normalized' );
subplot(2,1,1);stem(0:N,[A0*2,A(1:end)]); title('单边幅度谱');grid on;
subplot(2,1,2);stem(0:N,[Phi0,Phi(1:end)]);title('单边相位谱');grid on;

uf=fft(un,numel(un))/(Fs*T);    %因为每个周期有Fs个点
figure(4);
h=figure(4);
set(h,'position', [0.2 0.5 0.6 0.4],'Units','normalized' );  %以相对尺寸的方式绘图
set(h,'position', [0.2 0.5 0.6 0.4],'Units','normalized' );
set(h,'name','7.2.2 比较非周期信号傅立叶变换与周期信号傅立叶系数的关系');
subplot(2,1,1);stem(0:N,abs(uf(1:N+1)));title('非周期信号幅度谱（正频率部分）');grid on;
subplot(2,1,2);stem(0:N,angle(uf(1:N+1)));title('非周期信号相位谱（正频率部分）');grid on;
pause;

%7.2.3 用MATLAB符号运算函数fourier求解下列信号的傅立叶变换，并用ifourier验证

syms t;
ft1=exp(-t)*heaviside(t);%设置函数式
ft2=1/(sqrt(2*pi))*exp(-((t^2)/2));%设置函数式
Fw1=fourier(ft1);%对函数进行傅里叶变换
ft1i=ifourier(Fw1);%对函数进行傅里叶反变换
Fw2=fourier(ft2);
ft2i=ifourier(Fw2);
figure(7);
h=figure(7);
set(h,'position', [0.2 0.05 0.6 0.4],'Units','normalized' );
set(h,'position', [0.2 0.05 0.6 0.4],'Units','normalized' );
set(h,'name','7.2.3 用MATLAB符号运算函数fourier求解下列信号的傅立叶变换，并用ifourier验证');
subplot(2,1,1);ezplot(abs(Fw1));title('原信号1的幅度谱');axis([-6 6 0 1.1]); grid on;
subplot(2,1,2);h1=ezplot(abs(Fw2),[-6 6]);title('原信号2的幅度谱');axis([-6 6 0 1.1]); grid on;
set(h1,'Color','r','LineWidth',2)  %设置描述线为红色
figure(5);
h=figure(5);
set(h,'position', [0 0.5 0.5 0.4],'Units','normalized' );
set(h,'position', [0 0.5 0.5 0.4],'Units','normalized' );
set(h,'name','7.2.3 用MATLAB符号运算函数fourier求解下列信号的傅立叶变换，并用ifourier验证');
subplot(2,1,1);ezplot(ft1); title('原时域信号');axis([-0.1 3 0 1.1]); grid on;
subplot(2,1,2);h1=ezplot(ft1i); title('逆变换后的时域信号');axis([-0.1 3 0 1.1]); grid on;%通过傅里叶反变换得到的图像验证
set(h1,'Color','c','LineWidth',2);%设置描述线为蓝色
figure(6);
h=figure(6);
set(h,'position', [0.5 0.5 0.5 0.4],'Units','normalized' );
set(h,'position', [0.5 0.5 0.5 0.4],'Units','normalized' );
set(h,'name','7.2.3 用MATLAB符号运算函数fourier求解下列信号的傅立叶变换，并用ifourier验证');
subplot(2,1,1);h1=ezplot(ft2);title('原时域信号');axis([-3 3 0 0.5]); grid on;
set(h1,'Color','r');
subplot(2,1,2);h1=ezplot(ft2i);title('逆变换后的时域信号');axis([-3 3 0 0.5]); grid on;
set(h1,'Color','c','LineWidth',2);%通过傅里叶反变换得到的图像验证
pause;

%7.2.4 验证傅立叶变换的尺度变换性质

dt=0.1;
t1=-4:dt:4;
ft1=((1+cos(t1))/2).*(heaviside(t1+pi)-heaviside(t1-pi));%定义函数式ft1
ft2=((1+cos(2*t1))/2).*(heaviside(t1+0.5*pi)-heaviside(t1-0.5*pi));%定义函数式ft2
N=500;
k=-N:N;
W=pi*k/(N*dt);
Fw1=dt*ft1*exp(-1i*t1'*W);%分别对ft1和ft2进行傅里叶变换
Fw2=dt*ft2*exp(-1i*t1'*W);
figure(8);
h=figure(8);
set(h,'position', [0 0.1 0.5 0.8],'Units','normalized' );
set(h,'position', [0 0.1 0.5 0.8],'Units','normalized' );
set(h,'name','7.2.4 验证傅立叶变换的尺度变换性质');
subplot(2,1,1);plot(t1,ft1,'LineWidth',2); title('信号时域波形'); axis([-4 4 -0.1 1.1]); grid on;%绘制信号时域波形
subplot(2,1,2);plot(t1,ft2,'r','LineWidth',2); title('尺度变换后的时域波形');axis([-4 4 -0.1 1.1]); grid on;%绘制尺度变换后的信号时域波形
figure(9);
h=figure(9);
set(h,'position', [0.5 0.1 0.5 0.8],'Units','normalized' );
set(h,'position', [0.5 0.1 0.5 0.8],'Units','normalized' );
set(h,'name','7.2.4 验证傅立叶变换的尺度变换性质');
subplot(2,1,1);plot(W,real(Fw1),'LineWidth',2);title('信号的复合频谱');axis([-8 8 -0.6 3.5]); grid on;%绘制信号的复合频谱
subplot(2,1,2);plot(W,real(Fw2),'r','LineWidth',2);title('尺度变换后信号的复合频谱');axis([-8 8 -0.6 3.5]); grid on;%绘制尺度变换后的信号复合频谱
pause;

%7.3.2 幅度失真对听觉和视觉的影响

T=3;
Fs=40000;
t=0:1/Fs:T;
s1=sin(600*2*pi*t);%分别设置s1，s2，s3
s2=sin(850*2*pi*t);
s3=3*sin(850*2*pi*t);
s1s2=s1+s2;%定义函数之和的形式
s1s3=s1+s3;
figure(11); 
h=figure(11);
set(h,'position', [0.01 0.05 0.6 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.6 0.85],'Units','normalized' );
set(h,'name','7.3.2 幅度失真对听觉和视觉的影响');%绘制s1、s2、s3、s1s2、s1s3的时域图
subplot(511);plot(t(1:600),s1(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1');
subplot(512);plot(t(1:600),s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s2');
subplot(513);plot(t(1:600),s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s3');
subplot(514);plot(t(1:600),s1s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s2');
subplot(515);plot(t(1:600),s1s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s3');
[R,C]=meshgrid(0:1/Fs:600/160000:1/Fs:600/160000);
s1i=sin(600*2*pi*R);%分别设置s1i，s2i，s3i
s2i=sin(850*2*pi*R);
s3i=3*sin(850*2*pi*R);
s1s2i=sin(600*2*pi*R)+sin(850*2*pi*R);
s1s3i=sin(600*2*pi*R)+3*sin(850*2*pi*R);
figure(13);
h=figure(13);
set(h,'position', [0.62 0.05 0.37 0.85],'Units','normalized' );
set(h,'position', [0.62 0.05 0.37 0.85],'Units','normalized' );
set(h,'name','7.3.2 幅度失真对听觉和视觉的影响');
subplot(5,1,1),imshow(s1i(1:30,:)/3,[0,1]);title('s1i')
subplot(5,1,2),imshow(s2i(1:30,:)/3,[0,1]);title('s2i') 
subplot(5,1,3),imshow(s3i(1:30,:),[]);title('s3i') 
subplot(5,1,4),imshow(s1s2i(1:30,:),[]);title('s1s2i'); 
subplot(5,1,5),imshow(s1s3i(1:30,:),[]);title('s1s3i'); 
pause;
figure(11);
subplot(514);plot(t(1:600),s1s2(1:600),'r');axis([0,600/Fs,-3.7,3.7]);title('s1s2');
figure(13); %以声音的形式将信号的方式展现
subplot(5,1,4),imshow(mat2gray(s1s2i(1:30,:))*64,spring);title('s1s2i'); 
sound(s1s2,Fs);
pause(4);
figure(11);
subplot(514);plot(t(1:600),s1s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s2');
figure(13);
subplot(5,1,4),imshow(s1s2i(1:30,:),[]);title('s1s2i'); 
figure(11);
subplot(515);plot(t(1:600),s1s3(1:600),'r');axis([0,600/Fs,-3.7,3.7]);title('s1s3');
figure(13);
subplot(5,1,5),imshow(mat2gray(s1s3i(1:30,:))*64,spring);title('s1s3i'); 
sound(s1s3,Fs);
pause(4);
figure(11);
subplot(515);plot(t(1:600),s1s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s3');
figure(13);
subplot(5,1,5),imshow(s1s3i(1:30,:),[]);title('s1s3i'); 
pause;

%7.3.3 相位失真对听觉和视觉的影响

T=3;
Fs=40000;
t=0:1/Fs:T;
s1=sin(600*2*pi*t);%分别设置s1，s2，s3以及函数之和
s2=sin(850*2*pi*t);
s3=sin(850*2*pi*t-0.7);
s1s2=s1+s2;
s1s3=s1+s3;
figure(14); 
h=figure(14);
set(h,'position', [0.01 0.05 0.6 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.6 0.85],'Units','normalized' );
set(h,'name','7.3.3 相位失真对听觉和视觉的影响');
subplot(511);plot(t(1:600),s1(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1');
subplot(512);plot(t(1:600),s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s2');
subplot(513);plot(t(1:600),s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s3');
subplot(514);plot(t(1:600),s1s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s2');
subplot(515);plot(t(1:600),s1s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s3');
[R,~]=meshgrid(0:1/Fs:600/160000:1/Fs:600/160000);
s1i=sin(600*2*pi*R);%分别设置s1i，s2i，s3i以及函数之和
s2i=sin(850*2*pi*R);
s3i=sin(850*2*pi*R-0.7);
s1s2i=sin(600*2*pi*R)+sin(850*2*pi*R);
s1s3i=sin(600*2*pi*R)+3*sin(850*2*pi*R-0.7);
figure(16);
h=figure(16);
set(h,'position', [0.62 0.05 0.37 0.85],'Units','normalized' );
set(h,'position', [0.62 0.05 0.37 0.85],'Units','normalized' );
set(h,'name','7.3.3 相位失真对听觉和视觉的影响');
subplot(5,1,1),imshow(s1i(1:30,:),[0,1]);title('s1i')%分别绘制图像
subplot(5,1,2),imshow(s2i(1:30,:),[0,1]);title('s2i') 
subplot(5,1,3),imshow(s3i(1:30,:),[]);title('s3i') 
subplot(5,1,4),imshow(s1s2i(1:30,:),[]);title('s1s2i'); 
subplot(5,1,5),imshow(s1s3i(1:30,:),[]);title('s1s3i'); 
pause;
figure(14);
subplot(514);plot(t(1:600),s1s2(1:600),'r');axis([0,600/Fs,-3.7,3.7]);title('s1s2');
figure(16);
subplot(5,1,4),imshow(mat2gray(s1s2i(1:30,:))*64,spring);title('s1s2i'); 
sound(s1s2,Fs);
pause(4);
figure(14);
subplot(514);plot(t(1:600),s1s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s2');
figure(16);
subplot(5,1,4),imshow(s1s2i(1:30,:),[]);title('s1s2i'); 
figure(14);
subplot(515);plot(t(1:600),s1s3(1:600),'r');axis([0,600/Fs,-3.7,3.7]);title('s1s3');
figure(16);
subplot(5,1,5),imshow(mat2gray(s1s3i(1:30,:))*64,spring);title('s1s3i'); 
sound(s1s3,Fs);
pause(4);
figure(14);
subplot(515);plot(t(1:600),s1s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s3');
figure(16);
subplot(5,1,5),imshow(s1s3i(1:30,:),[]);title('s1s3i'); 
pause;

%7.4.2 对音频信号时移，与原信号在时域上和频域上进行比较

Fs=8000;
[testsou1,Fs1]=audioread('aud1.mp3');%读入音频
testsou1=testsou1(:,1);%对音频进行切分
testsou1=resample(testsou1,Fs,Fs1);%改变其维度
tem=testsou1(120000:120000+5*8000-1);
A1=tem;
A2=[zeros(Fs,1);A1];%时移一秒
A=[A1;zeros(Fs,1)];
Af=fft(A,Fs*5);  %使用快速傅里叶卷积函数
A2f=fft(A2,Fs*5);
figure(18);
h=figure(18);
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.2 对音频信号时移，与原信号在时域上和频域上进行比较');
subplot(2,1,1), plot([1:round(Fs*5/2)],abs(Af(1:round(Fs*5/2)))); title('Af')% 绘制其图像
subplot(2,1,2), plot([1:round(Fs*5/2)],abs(A2f(1:round(Fs*5/2)))); title('A2f')
figure(17);
h=figure(17);
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.2 对音频信号时移，与原信号在时域上和频域上进行比较');
subplot(2,1,1), plot(A), axis([0,numel(A), -1, 1]);grid on;title('A')% 绘制其图像
subplot(2,1,2), plot(A2), axis([0,numel(A2), -1, 1]);grid on;title('A2')
pause;
sound([A,A2],Fs);
subplot(2,1,1), plot(A,'r'), axis([0,numel(A), -1, 1]);grid on;title('A')% 绘制其图像
subplot(2,1,2), plot(A2,'r'), axis([0,numel(A2), -1, 1]);grid on;title('A2')
pause(6.6);
subplot(2,1,1), plot(A), axis([0,numel(A), -1, 1]);grid on;title('A')% 绘制其图像
subplot(2,1,2), plot(A2), axis([0,numel(A2), -1, 1]);grid on;title('A2')
pause;

%7.4.3 对音频信号进行尺度变换，与原信号在时域上和频域上进行比较

Fs=8000;
[testsou1,Fs1]=audioread('aud1.mp3');%读取音频
testsou1=testsou1(:,1);%对音频进行切分
testsou1=resample(testsou1,Fs,Fs1);%改变按采样率，实现抽取和内插的操作
tem=testsou1(120000:120000+5*8000-1);
A1=tem;
A3=resample(A1,numel(A1)*0.6,numel(A1));%尺度变换a=1/0.6
A1f=fft(A1,Fs*5); %使用快速傅里叶变换函数
A3f=fft(A3,Fs*5);
figure(20);
h=figure(20);
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.3 对音频信号进行尺度变换，与原信号在时域上和频域上进行比较');
subplot(2,1,1), plot([1:round(Fs*5/2)],abs(A1f(1:round(Fs*5/2)))); title('A1f')%对得到的图像进行绘制
subplot(2,1,2), plot([1:round(Fs*5/2)],abs(A3f(1:round(Fs*5/2)))); title('A3f')
figure(19);
h=figure(19);
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.3 对音频信号进行尺度变换，与原信号在时域上和频域上进行比较');
subplot(2,1,1), plot(A1), axis([0,numel(A1), -1, 1]);grid on;title('A1') %对得到的图像进行绘制

subplot(2,1,2), plot(A3), axis([0,numel(A1), -1, 1]);grid on;title('A3')
pause;
sound(A1,Fs);
subplot(2,1,1), plot(A1,'r'), axis([0,numel(A1), -1, 1]);grid on;title('A1')
pause(5)
subplot(2,1,1), plot(A1), axis([0,numel(A1), -1, 1]);grid on;title('A1') %对得到的图像进行绘制

pause(1.5)
sound(A3,Fs);
subplot(2,1,2), plot(A3,'r'), axis([0,numel(A1), -1, 1]);grid on;title('A3')
pause(4)
subplot(2,1,2), plot(A3), axis([0,numel(A1), -1, 1]);grid on;title('A3')
pause;

%7.4.4 对音频信号分别进行一次微分和二次微分操作，比较这两个信号与原信号的幅度谱和声音的变化。

Fs=8000;
[testsou1,Fs1]=audioread('aud1.mp3');%读入音频
testsou1=testsou1(:,1);%对音频进行切分
testsou1=resample(testsou1,Fs,Fs1);%进行抽样
tem=testsou1(120000:120000+5*8000-1);
A1=tem;
A4=[0;A1(1:end-1)]-A1; %求一次导数
A5=[0;A4(1:end-1)]-A4; %求二次导数
A1f=fft(A1,Fs*5);
A4f=fft(A4,Fs*5);
A5f=fft(A5,Fs*5);
figure(21);
h=figure(21);
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'name','对音频信号分别进行一次微分和二次微分操作');
subplot(3,1,1), plot([1:round(Fs*5/2)],abs(A1f(1:round(Fs*5/2)))); title('A1f')%对音频信号分别进行一次微分和二次微分操作并绘制图像
subplot(3,1,2), plot([1:round(Fs*5/2)],abs(A4f(1:round(Fs*5/2)))); title('A4f')
subplot(3,1,3), plot([1:round(Fs*5/2)],abs(A5f(1:round(Fs*5/2)))); title('A5f')
pause;
sound(A1,Fs);
subplot(3,1,1), plot([1:round(Fs*5/2)],abs(A1f(1:round(Fs*5/2))),'r'); title('A1f')
pause(6)
subplot(3,1,1), plot([1:round(Fs*5/2)],abs(A1f(1:round(Fs*5/2)))); title('A1f')
pause(0.5)
sound(A4,Fs);
subplot(3,1,2), plot([1:round(Fs*5/2)],abs(A4f(1:round(Fs*5/2))),'r'); title('A4f')
pause(6)
subplot(3,1,2), plot([1:round(Fs*5/2)],abs(A4f(1:round(Fs*5/2)))); title('A4f')
pause(0.5)
sound(A5,Fs);
subplot(3,1,3), plot([1:round(Fs*5/2)],abs(A5f(1:round(Fs*5/2))),'r'); title('A5f')
pause(6)
subplot(3,1,3), plot([1:round(Fs*5/2)],abs(A5f(1:round(Fs*5/2)))); title('A5f')
pause;

%7.4.5 对音频信号的傅立叶变换进行频移1Hz操作，比较其时域波形和声音的变化

Fs=8000;
[testsou1,Fs1]=audioread('aud1.mp3');%读入音频
testsou1=testsou1(:,1);%对音频进行切分
testsou1=resample(testsou1,Fs,Fs1);%取样
tem=testsou1(120000:120000+5*8000-1);
A1=tem;
Af=fft(A1);%使用快速傅里叶变换函数
Afsh=[zeros(round(numel(Af)/Fs),1);Af(1:round(numel(Af)/2));...
    Af(round(numel(Af)/2)+1:end);zeros(round(numel(Af)/Fs),1)];%频域1Hz
Ash=ifft(Afsh,40000);%频移后的时域信号，逆变换
figure(22);
h=figure(22);
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.5 对音频信号的傅立叶变换进行频移1Hz操作');
subplot(2,1,1), plot(A1), axis([0,numel(A1), -1, 1]);grid on;title('A1')
subplot(2,1,2), plot(real(Ash)), axis([0,numel(A1), -1, 1]);grid on;title('Ash')
figure(23);
h=figure(23);
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
plot([200:450],abs(Af(200:450)),'k');
hold on
plot([200:450],abs(Afsh(200:450)),'m','LineWidth',2);
grid on;
title('Af and Afsh')
pause;
sound(A1,Fs);
figure(22);
subplot(2,1,1), plot(A1,'r'), axis([0,numel(A1), -1, 1]);grid on;title('A1')%绘制其图像
pause(5)
subplot(2,1,1), plot(A1), axis([0,numel(A1), -1, 1]);grid on;title('A1')
pause(1)
sound(real(Ash),Fs);
subplot(2,1,2), plot(real(Ash),'r'), axis([0,numel(A1), -1, 1]);grid on;title('Ash')
pause(5)
subplot(2,1,2), plot(real(Ash)), axis([0,numel(A1), -1, 1]);grid on;title('Ash')
pause;

%7.5 看见声音

[testsou,Fs1]=audioread('aud1.mp3');%读入音频文件
testsou1=testsou(:,1);%对音频进行切分
Ig1=reshape(testsou1(1:100000),[250,400]); %对生成的矩阵重组为新矩阵
Ig1f=fft(testsou1,10000);%对其进行快速傅里叶变换
noise1=randn(100000,1);%生成噪声
ng1=reshape(noise1,[250,400]);%对噪声进行重组为新矩阵
ng1f=fft(noise1,10000);%对其进行快速傅里叶变换
figure(24);
h=figure(24);
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );%绘制音频灰度值图
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'name','7.5 看见声音');
subplot(2,2,1),imshow(Ig1,[]);title('Ig1')
subplot(2,2,2),imshow(ng1,[]);title('ng1')
subplot(2,2,3),plot(abs(Ig1f(1:5000)));title('Ig1f')
subplot(2,2,4),plot(abs(ng1f(1:5000)));title('ng1f')
pause;
sound([Ig1(:);zeros(40000,1);0.2*ng1(1:size(ng1,2)*size(ng1,1)*0.75)'],Fs1);  %顺序播放乐音和噪声
subplot(2,2,1),imshow(mat2gray(Ig1)*64,summer);title('Ig1');
pause(5);
subplot(2,2,1),imshow(Ig1,[]);title('Ig1');
pause(1);
sound(0.2*ng1(1:size(ng1,2)*size(ng1,1)*0.75),Fs1);  %顺序播放乐音和噪声
subplot(2,2,2),imshow(mat2gray(ng1)*64,summer);title('ng1');
pause(4);
subplot(2,2,2),imshow(ng1,[]);title('ng1');
pause;

%7.6 听见图像

function g=gscale(f,varargin)  %对gscale函数进行定义
if length(varargin)==0
    method='full8';
else method=varargin{1};
end
if strcmp(class(f),'double') & (max(f(:))>1 || min(f(:))<0)
    f=mat2gray(f);
end
switch method
    case 'full8'
        g=im2uint8(mat2gray(double(f)));
    case 'full16'
        g=im2uint16(mat2gray(double(f)));
    case 'minmax'
        low = varargin{2};
        high = varargin{3};
        if low>1 || low<0 || high>1 || high<0
            error('Parameters low and high must be in the range [0,1]')
        end
        if strcmp(class(f),'double')
            low_in=min(f(:));
            high_in=max(f(:));
        elseif strcmp(class(f),'uint8')
            low_in=double(min(f(:)))./255;
            high_in=double(max(f(:)))./255;
        elseif strcmp(class(f),'uint16')
            low_in=double(min(f(:)))./65535;
            high_in=double(max(f(:)))./65535;
        end
        g=imadjust(f,[low_in high_in],[low high]);
    otherwise
        error('Unknown method')
end

Fs=8000;
I1=imread('DSC_0333.JPG');%读入图片
%I2=imread('纹理7.JPG');
[R, C]=meshgrid(-2:0.01:2,-2:0.01:2);
u0=5*pi; v0=2*pi; A=1;
f1 = A*sin(u0*R + v0*C);%分别对f1，f2，f3进行定义
f2 = A*sin(u0*R) + A*sin(v0*C);
f3 = A*sin(u0*sqrt(R.^2+C.^2));
Ig2=gscale(f2);  %利用上面定义的函数对输入图像标度
Ig1=rgb2gray(I1);%将RGB图像转为灰度图像
Ig1=mat2gray(Ig1);%对图像矩阵的归一化
Ig1=imresize(Ig1,[400,800]);%修改尺寸
%Ig2=rgb2gray(I2);
Ig2=mat2gray(Ig2);
Ig2=imresize(Ig2,[400,800]);
Ig2fil=imfilter(Ig2,ones(15)/225); %采用滤波操作
Ig3=randn(400,800);
figure(25),
h=figure(25);
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'name','7.6 听见图像');
subplot(3,1,1),imshow(Ig1,[]);title('Ig1')
subplot(3,1,2),imshow(Ig2,[]);title('Ig2')
subplot(3,1,3),imshow(Ig3,[]);title('Ig3')
pause;
sound([Ig1(:);zeros(40000,1);Ig2fil(:);zeros(40000,1);0.2*Ig3(:)],Fs*8);
subplot(3,1,1),imshow(mat2gray(Ig1)*64,summer);title('Ig1')%绘制其图像
pause(5.6);
subplot(3,1,1),imshow(Ig1,[]);title('Ig1')
subplot(3,1,2),imshow(mat2gray(Ig2)*64,summer);title('Ig2')
pause(5.4);
subplot(3,1,2),imshow(Ig2,[]);title('Ig2')
subplot(3,1,3),imshow(mat2gray(Ig3)*64,summer);title('Ig3')
pause(5.5);
subplot(3,1,3),imshow(Ig3,[]);title('Ig3')

【实验感悟】

通过这次实验，我学会了对函数的对数幅频特性与相频特性绘制并对信号抽样与对抽样信号的频谱分析，以及对抽样后的信号进行重建。并且通过实验，我加深了对于信号处理这门课的理解，通过形象的方式了解到信号的特点，其中有音频信号，也有图像信号。对于信号的重构让我看到了乐音与噪声的不同之处，通过实验让我更加形象地理解了任何周期信号都可以分解成正弦信号加权的形式。而且本次实验也让我从不同的维度上（音频、图像等等）对于高频信号与低频信号有了更加深刻的理解。此外，我还通过查找文档文献的方法掌握了与信号处理相关的函数的使用方法，例如 imshow(I,[low high])函数用来显示灰度图像I,以二元素向量 [low high] 形式指定显示范围，y = resample(x,p,q)使用多相滤波器实现对矢量X中的序列在原始采样率的P/Q倍上重新采样，reshape函数将原矩阵重组为新矩阵，fft快速傅里叶卷积函数等等，这也让我从工程技能上有所收获。

本文参与腾讯云自媒体同步曝光计划，分享自作者个人站点/博客。

原始发表：2022-07-04，如有侵权请联系 cloudcommunity@tencent.com 删除

matlab