matlab信号处理学习（实战代码） |
您所在的位置：网站首页 › sq在电路 › matlab信号处理学习（实战代码） |
matlab信号处理学习（实战代码）
|
没有下载matlab可以打开网页版Octave，很好用Octave Online · Cloud IDE compatible with MATLAB (octave-online.net) part1 %创建正弦波 %定义信号采样序列。从0s到1s每隔0.001s采样一次，共采样1000次 t=0:0.001:1; y1=sin(2*pi*50*t); y2=sin(2*pi*50*t)+2*sin(2*pi*100*t); y3=y2+0.5*randn(size(t));%加一个噪声信号 %这里的两个50要保持一致，50是取50个采样，最多取t中定义的1000个 subplot(131),plot(t(1:50),y1(1:50)),title('50HZ正弦波') subplot(132),plot(t(1:50),y2(1:50)),title('50HZ+100HZ正弦波') subplot(133),plot(t(1:20),y3(1:20)),title('50HZ+100HZ噪声正弦波')%创建周期性方波 t=linspace(0,3*pi)'; x=square(t);%利用函数创建方波信号x plot(t,x); xlabel('t/\pi')%t/π,横轴 xlabel('t/pi')%t/pi %方波是一种非正弦曲线的波形，震荡电路出来的正弦波一般都有谐波分量，方波就是由一系列的谐波分量叠加而成，理想的方波只有0-1两个值%创建非周期性方波 t=-4:0.001:4; y=rectpuls(t);%产生非周期方波信号，默认方波宽度为1 z=square(t);%产生周期性的矩形脉冲信号 subplot(131) plot(t,z,'b','linewidth',4); axis([-4,4,-1.5,1.5]); grid on; xlabel t,ylabel n(t); subplot(132) plot(t,y,'r','linewidth',4) axis([-4,4,-1.5,1.5]); grid on; xlabel t,ylabel h(t); y=2*rectpuls(t,2);%产生指定宽度为2的非周期方波信号。前面的*2是高度*2 subplot(133) plot(t,y,'linewidth',4) axis([-4,4,-0.5,2.5]); grid on;%创建锯齿波 f=50; T=5*(1/50);%采样时间0.1s，信号周期为5，频率为50HZ fs=1000;%定义采样率为1000 t=0:1/fs:T-1/fs; y=sawtooth(2*pi*50*t,0);%周期为2pi，0是向左倾斜的锯齿波 z=sawtooth(2*pi*50*t,1);%周期为2pi，1是向右倾斜的锯齿波 subplot(211); plot(t,y);%创建指数波 fs=100; t=0:1/fs:1; C=1; a=1; y=C*exp(a*t); subplot(211); plot(t,y); C=-1; a=-1; y=C*exp(a*t); subplot(212); plot(t,y,'.');%创建斜坡信号,线性增长信号 fs=10; t=0:1/fs:1; x=t; plot(t,x,'r'); hold on stem(t,x,'bo');%增状图命名stem xlabel('时间序列')； ylabel('x(t)'); title('斜坡序列')%创建插值序号——解决由于信号采样频率降低发生的频率混叠 fs=100; t=0:1/fs:2-1/fs; x=sin(2*pi*10*t.^3).^2;%外侧.^2变密,内侧.^3各个幅度变乱不规则 y=interp(x,4);%输入插值信号函数,将插值信号采样率增加一个因子4 tiledlayout(2,1); nexttile; stem(t,x,'filled',markersize,3) grid on xlabel('Sample Number') ylabel('Original') nexttile;%跟subplot的作用一样 t=1:length(y); stem(t,y,'filled','markersize',1) xlabel('Sample Number') ylabel('Interpolated')%创建添加噪声的正弦波 fs=1000; t=0:1/fs:1-1/fs; y1=sin(2*pi*t); y2=wgn(1000,1,0);%创建1000*1的高斯白噪声y2,噪声样本功率为0 y3=sin(2*pi*t)+wgn(1000,1,0);%创建添加噪声的正弦波 y4=y1+0.5*randn(size(t));%在正弦波上添加正态分布的随机噪声信号，噪声信号的幅值是正弦波的一半 subplot(221),plot(t,y1),title('正弦波信号'); subplot(222),plot(t,y2),title('高斯白噪声信号'); subplot(223),plot(t,y3),title('添加高斯白噪声的正弦波信号'); subplot(224),plot(t,y4),title('添加正态分布噪声的正弦波信号');%创建添加噪声的指数波 fs=1000; t=0:1/fs:1-1/fs; x=airy(t*10).*exp(-t,^2); y1=x+wgn(1000,1,1);%添加高斯白噪声，噪声功率为1 y2=awgn(x,1,'measured');%高斯白噪声信噪比为1，measured指在加入噪声前测定信号强度 subplot(311),plot(t,x),title('指数波信号'); subplot(312),plot(t,y1),title('叠加高斯白噪声的指数波信号'); subplot(313),plot(t,y2),title('添加高斯白噪声的指数波信号');%创建随机波 t=0:100; N=length(t); x=rand(1,N);%创建随机信号 y=randn(1,N);%创建正态分布的随机信号 subplot(131); plot(t,x,'k'); ylabel('x(t)'); subplot(132); plot(t,y,'r'); ylabel('y(t)'); subplot(133); stem(t,x,'filled','g'); ylabel('x2(t)');%创建sinc波，sinc(t)=sin(t)/t x=linspace(-5,5); y1=sin(pi*x)./(pi*x); y2=sinc(x); plot(x,y1,'o',x,y2); xlabel Time,ylabel Signal; %在图形中添加图例便于区分曲线 legend('sin signal','sinc signal','Location','NorthWest') legend boxoff;%图例不显示边框%创建对数扫频信号 t=0:0.001:1; y=chirp(t,1e-6,1,50,'logarithmic');%初始时刻瞬时频率1e-6，t=1时，瞬时频率为50HZ，扫频方法为logarithmic的对数扫频 plot(t,y); axis([0.8,1,-1,1]);%调整坐标范围%创建迪利克雷信号diric x=0:0.001:4*pi; subplot(211); plot(x,diric(x,7));%计算x的7级迪利克雷函数，并绘制xinhaobox axis tight,title('n=7'); subplot(212); plot(x/pi,diric(x,8)); title('n=8'),xlabel('x/\pi');%创建三角波信号 t=-3:0.1:3; %采样时间6s y=tripuls(4*t,2*pi,0.2);%采样次数为4t，周期为2pi,斜率为0.2的非周期三角波 z=sawtooth(4*t,0);%周期为2pi的周期性三角波（锯齿波） subplot(121); plot(t,y,'r^','linewidth',2); grid on; subplot(122); plot(t,z,'linewidth',2); grid on;part2 啦啦啦啦部分加上跑的图啦 %熟悉二维图形的绘制的一些工具 %显示4*4图形分割 t1=(0:11)/11*pi; t2=(0:400)/400*pi;%数字越大，越密集 t3=(0:50)/50*pi; y1=cos(t1).*cos(5*t1); y2=cos(t2).*cos(5*t2); y3=cos(t3).*cos(5*t3); subplot(221),plot(t1,y1,'r.'); title('(1)点过少的离散图形'); subplot(222),plot(t1,y1,t1,y1,'r.'); title('(2)点过少的连续图形'); subplot(223),plot(t2,y2,'r.'); title('(3)点过多的离散图形'); subplot(224),plot(t3,y3); title('(4)实线画连续图形');%图窗布局的应用 x=linspace(-pi,pi); y=sin(x); tiledlayout(2,2);%2*2共四个图块 nexttile;%激活第一个图块 plot(x); nexttile;%激活第一个图块 plot(x,y); nexttile([1,2]);%创建第三个图块，该图块占据一行两列的坐标区 plot(x,y);%保持命令的应用 x=linspace(-pi,pi); y1=sin(x).*exp(x); plot(x,y1); hold on;%使得第二条曲线叠加在第一条曲线上面，而不是直接覆盖 y2=cos(x).*sin(x); plot(x,y2); hold off %关闭图像保持命令，此时绘制第三条曲线就会覆盖掉前两个 y3=2*sin(3*x); plot(x,y3)%调整坐标系 x=0:pi/100:2*pi; y=cos(x).^5; plot(x,y); axis([0,pi,-2,2]);%绘制三角函数 x=linspace(0,10*pi,100); plot(x,cos(x).^2+sin(x)); title('三角函数之和'); xlabel x坐标,ylabel y坐标%绘制函数图形并标注 x=linspace(-5,5); y=x.^3+exp(x); plot(x,y); xt=[-3,3];%设置标注位置 yt=[-24,46];%设置标注位置 str={'loccal min','local max'}; text(xt,yt,str);%添加标注%绘制向量图形，并标出函数名 plot(1:10); gtext('My Plot','Color','red','FontSize',14);%添加绘图注释 x=0:0.1:5; y1=exp(0.5*x).*cos(2*x); y2=cos(x).^5; plot(x,y1,'r-',x,y2,'b-.'); title('函数曲线'); %添加名字 legend('函数1','函数2'); %添加图例 xlabel 自变量x,ylabel 函数值y; grid on;%显示网格part3 %傅里叶变换与谱分析！！！本节有很多matlab'中的函数没有在octave中实现！以下代码建议在matlab中跑！ 按照自己的理解，大部分逆变换都相当于对变换的撤销动作  %信号变换的本质上将信号从时域转化为频域
fs=100;
N=1000;
n=0:N-1;
t=n/fs;
x=diric(t,7);%随时间变化的迪利克雷信号
subplot(311);
plot(t,x);
title('原始函数');
subplot(312);
y=fft(x);%快速傅里叶变换
plot(t,y);
title('傅里叶变换');
subplot(313);
z=ifft(y);
plot(t,z);
title('傅里叶逆变换');%发现z函数和x函数一样 %计算周期方波信号与非周期方波信号的频域变换
fs=100;
t=0:1/fs:2;
x=[square(2*pi*t);rectpuls(2*pi*t)];%定义函数矩阵,两个元素分别是周期方波信号与非周期方波信号
subplot(311);
plot(t,x,'linewidth',4);
axis([0 2 -1.5 1.5]);
title('原始信号');
subplot(312);
y=fft2(x);%二维快速傅里叶变换
plot(t,y);
title('傅里叶变换后信号');
axis([0 2 -11 20]);
subplot(313);
z=ifft2(y);
plot(t,z);
title('二维傅里叶逆变换');%一维傅里叶变换是正弦信号的叠加，而二维傅里叶变换是正弦平面波的叠加。
axis([0 2 -1.5 1.5]); 在matlab中，经过fft变换后，数据的频率范围是从[0,fs]排列的。而一般，我们在画图或者讨论的时候，是从[-fs/2,fs/2]的范围进行分析。因此，需要将经过fft变换后的图像的[fs/2,fs]部分移动到[-fs/2,0]这个范围内。而fftshift就是完成这个功能。通常，如果想得到所见的中间是0频的图像，经过fft变换后，都要再经过fftshift。  %计算余弦波噪声信号的傅里叶变换
fs=100;%信号采样率100HZ
t=0:1/fs:1-1/fs;%信号采样时间序列
f=1;%设置信号频率为1
x=cos(2*pi*f*t)+rand(1,100);%在余弦信号中添加随机噪声
subplot(311);
plot(t,x);
title('原始信号');
subplot(312);
y=fftshift(x);%对信号x进行0频分量转移计算
plot(t,y);
title('傅里叶变换')
subplot(313);
z=ifftshift(y);%对信号x进行0频分量转移逆运算
plot(t,z);
title('傅里叶逆变换'); %(N维傅里叶变换)对信号进行fft变换
fs=100;%信号采样率100HZ
f1=100;f2=200;%定义两个信号的频率
t=0:1/fs:1;%信号采样时间序列
x1=sin(2*pi*100*t)+cos(2*pi*200*t);%原始x
subplot(311);
plot(x1);
title('原始信号');
x=sin(2*pi*100*t)+cos(2*pi*200*t)+rand(size(t));%x收到噪声污染
subplot(312),plot(x(1:50));%绘制时域内的信号，噪声污染后很难看出正弦波的成分
title('噪声污染信号');
Y=fft(x,512);%对x进行512点的傅里叶变换
f=fs*(0:256)/512;
subplot(313),plot(f,Y(1:257));%绘制频域内信号，可以明显看出信号中100HZ和200HZ的两个频率分量
title('加噪声信号频域变换'); %离散傅里叶变换 N=20;%定义信号采样点数 n=0:N-1;%定义信号采样时间序列 x=exp(n);%定义原始指数信号 D=dftmtx(N);%求DFT矩阵 X=x*D; subplot(211); stem(n,x),title('原始信号x'); subplot(212); stem(abs(X));title('离散傅里叶变换X'); %复倒谱分析 matlab跑
%复倒谱是指一个信号的傅里叶变换对数的傅里叶反变换
fs=1000;
t=0:1/fs:1-1/fs;
x=square(20*pi*t);%创建方波信号
x=vco(x,[0.1 0.5]*fs,fs);%对方波信号进行调频，第二个参数是调频范围！！！！！！！！注意此处VCO功能属于Octave Forge的信号包，但尚未实现
y=fft(x);%快速傅里叶变换
z=cceps(x);%复倒谱分析
ax1=subplot(311),plot(t,x);
title('原始信号');
ax2=subplot(312);
plot(t,y);
title('原始信号进行fft');
ax3=subplot(313);
plot(t,z);
title('原始信号复倒谱');%傅里叶同步压缩变换 在matlab中跑(m跑)
fs=1000;
t=0:1/fs:1;
x=square(20*pi*t)+0.2*randn(size(t));%受污染的方波信号
subplot(211),plot(t,x);
title('原始信号');
subplot(212);
fsst(x,fs,'yaxis');%计算并绘制信号的傅里叶同步压缩变换
title('原始信号傅里叶同步压缩clc变换');%对信号进行傅里叶同步压缩变换和逆变换 m跑
fs=1000;
t=0:1/fs:3;
x=chirp(t,30,2,5).*exp(-(2*t-3).^2)+2;%创建高斯调制啁啾信号，信号直流值为2，初始啁啾频率为30HZ，2s后衰减到5Hz
subplot(311),plot(t,x);%绘制时域内原始信号
xlabel('Time(s)');
title('原始信号x');
subplot(312);
[y,f]=fsst(x,fs); %计算信号的傅里叶同步压缩变换
fsst(x,fs); %绘制信号的傅里叶同步压缩变换
xlabel('Time(s)');
title('原始信号傅里叶同步压缩变换');
z=ifsst(y);%进行傅里叶同步压缩逆变换，反转变换重建信号
t=(0:length(y)-1)/fs; %定义信号变换后的采样时间序列
subplot(313),plot(t,z);%绘制重构信号
xlabel('Time(s)');
title('重构信号');
%结果是第一个和第三个信号一致%计算信号dft变换（使用goertzel算法）
%goertzel算法目的是从给定的采样中求出某一特定频率信号的能量，用于有效性评价
fs = 1000;
f1=10e3;f2=1.24e3;
t=0:1/fs:1;
x= cos(2*pi*t*f1)+cos(2*pi*t*f2)+randn(size(t));%受污染的余弦波信号
subplot(2,1,1),plot(t,x);
title('原始信号x')
N = (length(x)+1)/2; %计算采样点数
f = (fs/2)/N*(0:N-1); %计算频率
indxs = find(f>1.2e3 & f 瀑布图%对信号执行短时傅里叶变换和逆变换
fs = 1000;
t = 0:1/fs:2-1/fs;
f0 = 100;%初始瞬时频率
f1 = 500;%参考时间瞬时频率
x =chirp(t,f0,1,f1,'quadratic',[],'concave');%先行扫频余弦信号，扫频方法为quadratic二次扫频，[]表示忽略初始相位，谱图形状为凸抛物线
win = hamming(100,'periodic'); %设计一个长度为100的周期性hamming窗
noverlap = 80; %样本重叠数为80
subplot(3,1,1)
plot(t,x)
xlabel('时间t');ylabel('原始信号x(t)');
title('原始信号');
subplot(3,1,2)
[y,f,t] = stft(x,fs,'Window', win,'OverlapLength',noverlap);%使用重叠样本数80计算信号的短时傅里叶变换
plot(f,y) ;
xlabel('频率f');ylabel('变换信号y(t)');
title('STFT变换');
subplot(3,1,3)
[iy,t] = istft(y,fs,'Window', win,'OverlapLength',noverlap);
plot(t,iy) %使用重叠样本数80计算信号的短时傅里叶逆变换
xlabel('时间t');ylabel('变换信号iy(t)');
title('ISTFT变换');%绘制短时傅里叶变换谱图spectrogram
t = -2:1/1e3:2;%利用线性间隔值函数定义信号采样时间序列，信号持续时间4s
f = 10;
y = sin(2*pi*f*t.^3);
spectrogram(y);%对比 窗函数与长度 对信号的短时傅里叶变换谱图 的影响
%加窗实质是用一个所谓的窗函数与原始的时域信号作乘积的过程（当然加窗也可以在频域进行，但时域更为普遍），使得相乘后的信号似乎更好地满足傅立叶变换的周期性要求。
fs = 3000;
t = 0:1/fs: 1-1/fs;
x= chirp(t,100,1,400,'quadratic') + chirp(t,350,1,50);%创建叠加的先行扫频余弦信号
%创建四个长度不同的taylorwin窗口
window1= taylorwin(256);
window2= taylorwin(512);
window3= taylorwin(1024);
window4= taylorwin(2056);
%创建四个长度不同的矩形窗口
window11= rectwin(256);
window21= rectwin(512);
window31= rectwin(1024);
window41= rectwin(2056);
tiledlayout(5,2)%五行两列的分布图布局
nexttile;
spectrogram(x,32);%将信号分为32个样本段，计算并绘制信号的短时傅里叶变换谱图
nexttile;
spectrogram(x,window1,'yaxis');%利用长度为256泰勒窗计算并绘制信号的短时傅里叶变换谱图，指定频率轴绘制在y轴
nexttile;
spectrogram(x,window1);
nexttile;
spectrogram(x,window2);
nexttile;
spectrogram(x,window3);
nexttile;
spectrogram(x,window4);
nexttile;
spectrogram(x,window11);
nexttile;
spectrogram(x,window21);
nexttile;
spectrogram(x,window31);
nexttile;
spectrogram(x,window41);
%结论：对于短时傅里叶变换，最重要的是窗口长度的选取。当频域刻度和平移步长足够密时，增加的只是生成图像的大小，但是物理层面的分辨率没有改变。改变物理层面分辨率的是窗口长度
%窗口长度大：频率清晰，时间模糊
%窗口长度小：频率模糊，时间清晰 %对比重叠样本数对短时傅里叶变换谱图的影响
fs = 1000;
t = 0:1/fs: 1-1/fs;
f = 2;
x= exp(1j*pi*sin(2*pi*f*t)*32);
noverlap1=10; %定义不同的信号重叠样本数
noverlap2=20;
tiledlayout(3,1)
nexttile;
spectrogram(x,32,noverlap1,64,'centered','yaxis');%分成32个样本段绘制信号的短时傅里叶变换中心双边谱图，重叠样本数10（重叠样本数必须小于样本段数32，频率显示在y轴）
nexttile;
spectrogram(x,32,noverlap2,64,'centered','yaxis'); %重叠样本数20，对比重叠样本数对频谱的影响
nexttile;
spectrogram(x,32,noverlap1,64,'centered','reassigned','yaxis');%计算并绘制信号的短时傅里叶变换中心双边谱图，将能量集中 %绘制短时傅里叶变换单边与双边谱图
fs = 1000;
t = 0:1/fs: 1-1/fs;
f = 2;
x= chirp(t,100,1,400,'quadratic');%扫频余弦信号
Nsample=32; %样本分割段数
noverlap=10; %信号重叠样本数
nfft=64; %DFT点数
tiledlayout(3,1)
nexttile;
spectrogram(x,Nsample,noverlap,64,'onesided','yaxis'); %onesided将信号分成32个样本段，绘制信号的短时傅里叶变换单边谱图，频率显示在y轴
nexttile;
spectrogram(x,Nsample,noverlap,64,'centered','yaxis');%centered将信号分成32个样本段，绘制信号的短时傅里叶变换中心双边谱图
nexttile;
spectrogram(x,Nsample,noverlap,64,'twosided','yaxis'); %twosided估计信号重新分配的双边谱图 %创建对称凸二次啁啾信号
fs=1000;
t = -2:1/fs:2;
f0 = 100;
f1 = 200;
x = chirp(t,f0,1,f1,'quadratic',[],'convex');
y = vco(x,[0.1 0.8]*fs,fs); %调频
tiledlayout(2,2)
nexttile;
pspectrum(x,t,'spectrogram','TimeResolution',0.1, ...
'OverlapPercent',99,'Leakage',0.85)%计算并绘制啁啾信号的谱图。指定相邻段之间99%的重叠和0.85的频谱泄露。
title('频谱图')
nexttile;
stft(x,fs); %计算并绘制信号的STFT
title('短时傅里叶变换')
nexttile;
spectrogram(x,100,80,1024,fs,'yaxis') %计算并绘制啁啾的短时傅里叶变换谱图。将信号分成100个样本段，指定相邻段之间重叠的80个样本、DFT长度为1024个样本
title('短时傅里叶变换谱图')
nexttile;
xspectrogram(x,y,100 ,80,1024,fs,'yaxis') %计算并绘制啁啾的交叉功率谱图。将信号分成100个样本段，指定相邻段之间重叠的80个样本、DFT长度为1024个样本
title('短时傅里叶变换交叉功率谱图') part4 频谱估计%采样点数对正弦波信号的傅里叶变换的影响
fs = 100;
N1=100;
N2=200;
n1 = 0:N1-1;
n2 = 0:N2-1;
t1=n1/fs;
t2=n2/fs;
x1 = sin(2*pi*4*t1);%两个正弦波信号
x2 = sin(2*pi*4*t2);
subplot(3,1,1)
plot(t1,x1,'*',t2,x2,'r')
xlabel('时间/s');
legend('N1=100','N2=200');
title('原始函数')
grid on;
subplot(3,1,2)
y1=fft(x1,N1);
y2=fft(x2,N2);
plot(t1,y1,t2,y2)
xlabel('时间/s');
title('傅里叶变换')
legend('N1=100','N2=200');
grid on;
subplot(3,1,3)
z1=abs(y1);%计算信号1傅里叶变换后的振幅
z2=abs(y2);%计算信号2傅里叶变换后的振幅
f1=(0:N1-1)*fs/N1;%定义频率序列
f2=(0:N2-1)*fs/N2;
plot(f1(1:N1/2),z1(1:N1/2)*2/N1,f2(1:N2/2),z2(1:N2/2)*2/N2);
xlabel('频率/Hz');ylabel('振幅');
title('傅里叶幅值变换')
legend('N1=100','N2=200');
grid on; %绘制正弦波形的单边和双边频谱
f0 = 100;
fs = 5000;
n=1:1:1000;
N = length(n);
t =n/fs;
x = sin(2*pi*f0*t);
subplot(3,1,1)
plot(n,x)
title('原始信号')
y =fft(x);
y2=abs(y/N); %计算信号双边谱的幅度！！！！公式记得
subplot(3,1,2)
plot(n,y2)
xlabel('f(Hz)')
ylabel('|P2(f)|')
title('FFT信号双边频谱')
y1=y2(1:N/2+1);
%实际上，双边谱是由单边谱/2所得频谱，然后正负对称频率叠加到一起，而零频率不需要/2
y1(2:end-1)=2*y1(2:end-1); %计算信号单边谱幅度！！！！公式记得
f = fs*(0:N/2)/N; %计算信号频谱幅度所对应的频率
subplot(3,1,3)
plot(f,y1)
xlabel('f(Hz)')
ylabel('|P1(f)|') %信号单边频谱
title('FFT信号单边频谱') %绘制正弦信号的相位谱
fs=1000; %定义采样序列为1000Hz
N=1024; %定义采样点数
n=0:N-1; %定义采样点数序列
t=n/fs; %定义采样时间序列
y=sin(2*pi*100*t); %定义正弦信号
Y=fft(y,N); %计算FFT信号
A=abs(Y); %求FFT信号的幅值
f=n*fs/N; %计算信号频率
subplot(311)
%奈奎斯特频率：为了防止信号混叠，需要定义的最小采样频率，此频率称为奈奎斯特频率。奈奎斯特频率是离散信号系统采样频率的一般f_nyquist=fs/2.根据对称性，用fft做信号频谱分析，只需要考察0~f_nyquist范围内的幅频特性即可f_nyquist=fs*(0:N/2)/N
plot(f(1:N/2),A(1:N/2)); %绘制奈奎斯特频率相位谱
xlabel('频率/hz'),ylabel('幅值'),title('幅值谱');
grid on;
subplot(312)
ph=2*angle(Y(1:N/2)); %求信号的相位
plot(f(1:N/2),ph(1:N/2)); %绘制频率相位谱
xlabel('频率/hz'),ylabel('相位'),title('相位谱');
grid on;
subplot(313)
ph1=ph*180/pi; %将信号相位由弧度转化为角度
plot(f(1:N/2),ph1(1:N/2)); %绘制频率相位谱
xlabel('频率/hz'),ylabel('相位'),title('相位谱');
grid on;
 %创建余弦波FFT相位谱
fs = 100;
t = 0:1/fs: 4-1/fs;
f1=10;f2=30;f3=40;
x=cos(2*pi*f1*t)+2*cos(2*pi*f2*t-pi/4)+4*cos(2*pi*f3*t+pi/2); %创建3个余弦波组成的信号。第一个余弦波相位是0，幅值是1；第二个余弦波相位为-pi/4，幅值是2；第三个余弦波相位为pi/2，幅值是4
x= awgn(x,10,'measured'); %创建添加高斯白噪声的余弦波，高斯白噪声功率值为10
subplot(221)
lx = length(x); %获取傅里叶变换后信号的长度
f = (-lx/2:lx/2-1)/lx*fs;
stem(f,abs(x)),title('原始信号序列图')%计算信号的幅值，绘制频率-幅值函数
xlabel 'Frequency (Hz)'
ylabel '|x|'
grid
y = fft(x);
subplot(222)
stem(f,abs(y)),title('FFT信号序列图')%计算FFT信号的幅值，绘制频率-幅值函数
xlabel 'Frequency (Hz)'
ylabel '|y|'
grid
subplot(223)
z = fftshift(y);%将傅里叶变换后信号中的零频分量移到频谱中心
stem(f,abs(z)),title('FFT频谱序列图') ;%计算信号的傅里叶变换，将变换幅值绘制为频率函数
xlabel 'Frequency (Hz)'
ylabel '|z|'
grid
subplot(224)
tol=1e-6; %定义小幅值变换值
z(abs(z)%光谱泄露对正弦信号频谱的影响
fs = 100;
t = 0:1/fs:2-1/fs;
x=sin(2*pi*20*t);
subplot(311)
pspectrum(x,t,'Leakage',1) %将泄漏增加到最大值1，解决了紧密间隔的音调，但掩盖了附近的微弱音调
subplot(312)
pspectrum(x,t,'Leakage',0) %将泄漏降低到最小值0，以牺牲光谱分辨率为代价，将泄漏降低到最小
subplot(313)
pspectrum(x,t,'Leakage',0.5) %将泄露显示默认值0.5 %信号功率谱单边谱与双边谱（三角波）
f = 50;
T = 10*(1/f);%10个信号周期
fs = 1000;
t = 0:1/fs:T-1/fs;
x = sawtooth(2*pi*f*t,1/2);%产生基频为50Hz的三角波的10个周期
subplot(221)
plot(t,x)
title('三角波信号')
subplot(222)
pspectrum(x,fs,'Leakage',0.91) %默认情况下绘制功率谱的单边谱，泄漏量0.91
subplot(223)
pspectrum(x,fs,'Leakage',1,'TwoSided',true) %TwoSided绘制双边光谱，泄漏量为1
subplot(224)
[p,f,t] = pspectrum(x,fs,'spectrogram'); %默认情况下绘制单边谱
waterfall(f,t,p'); %瀑布图形式显示光谱
xlabel('Frequency (Hz)')
ylabel('Time (seconds)')
view([60 45])
colormap hot %颜色设置为JEF %创建对称凸二次啁啾波频谱
fs=1e3;
t = -2:1/fs:2;
f0 = 100;
f1 = 200;
x = chirp(t,f0,1,f1,'quadratic',[],'convex');%余弦扫频信号
subplot(221)
pspectrum(x,t,'spectrogram','TimeResolution',0.1, ...
'OverlapPercent',99,'Leakage',0.85) %计算绘制啁啾谱图，将信号分割成分段，时间分辨率为0.1s，相邻段之间重叠为99%，频谱泄漏量为0.85
subplot(222)
pspectrum(x,fs,'spectrogram', ...
'FrequencyLimits',[100 300],'TimeResolution',1)%计算信号谱图，频带设置在100~300，谱图时间分辨率为1
subplot(223)
pspectrum(x,fs,'FrequencyLimits',[100 300],'Leakage',1) %计算信号功率谱，频带设置在100~300，频谱泄露量为1
subplot(224)
pspectrum(x,fs,'persistence', ...
'FrequencyLimits',[100 300],'TimeResolution',1) %计算信号的持续谱，频带设置在100~300，谱图时间分辨率为1 %绘制信号谱图
fs = 1000;
t=0:1/fs:1;
x= exp(2j*pi*100*cos(2*pi*2*t)) + randn(size(t))/100;%受噪声污染的指数波信号
subplot(2,1,1)
[sp,fp,tp] = pspectrum(x,fs,'spectrogram','FrequencyResolution',10); %指定频率分辨率为10Hz，计算光谱图[光谱，频谱频率，时间值]
mesh(tp,fp,sp) %三位网格图命令绘制光谱图
xlabel('Time (s)')
ylabel('Frequency (Hz)')
subplot(2,1,2)
waterfall(fp,tp,sp');%绘制光谱瀑布图
xlabel('Frequency (Hz)')
ylabel('Time (seconds)')
view([30 45])%调整视图视角 %创建信号的谱图和重分配谱图
fs = 1000;
t=0:1/fs:1;
x= (1+0.5*cos(2*pi*5*t)).*cos(2*pi*50*t+5*sin(2*pi*10*t));%余弦波信号
x=hilbert(x);%对信号进行希尔伯特变换，返回离散时间解析信号
subplot(2,2,1)
pspectrum(x,fs,'spectrogram'); %计算绘制光谱图
subplot(2,2,2)
pspectrum(x,fs,'spectrogram','FrequencyResolution',10,'Reassign',true) %指定频率分辨率带宽10Hz，计算绘制重新分配的光谱图
subplot(2,2,3)
pspectrum(x,fs,'spectrogram','TimeResolution',0.2) %用0.2s的时间分辨率计算光谱图
subplot(2,2,4)
pspectrum(x,fs,'spectrogram','TimeResolution',0.2,'Reassign',true)%用0.2s的时间分辨率计算绘制重新分配的光谱图 %功率谱几个概念： %1.信号以波的形式存在着，会产生功率。单位频带的信号功率称为功率谱，它显示在一定区域中信号功率随着频率变化的分布情况。 %2.对于确定性的信号，特别是非周期的确定性信号，常用能量谱描述。对随机信号，无法用确定的时间函数来表示，也就无法用频谱表示，往往用功率谱描述其频率特性 %3.功率谱是功率谱密度函数的简称，是单位频率的信号功率（单位频率内的信号能量）。表示了信号功率随着频率的变化情况； %能量谱几个概念 %1.对于任意的时间信号x(t),这个信号可以是任意随时间变化的物理量。在对信号进行能量分析时，不加区分的将其视为 1欧姆的电阻上的电流。这个单位电阻的能量属性，就是这个信号的能量属性 %2.如果信号是能量信号，通过傅里叶变换，就很容易分离不同频域分量所对应的能量。 %功率谱密度PSD的几个概念 %1.定义了信号或时间序列的功率如何随频率分布 %2.cpsd命令用于计算信号的交叉功率谱密度 %创建信号交叉功率谱 fs = 1000; t = 0:1/fs: 1-1/fs; f0 = 200; x = sin(2*pi*f0*t).^2; subplot(3,1,1) plot(t,x) %绘制时域图 xlabel('时间序列t');ylabel('原始信号x(t)'); title('原始信号'); subplot(3,1,2) pspectrum(x,t,'Leakage',1) %绘制泄漏增为最大值的品偶 title('信号频谱'); subplot(3,1,3) r = wgn(1000,1,0);%创建高斯白噪声 cpsd(x,r,triang(500),128,2048,fs); %绘制信号三角窗分割的交叉功率谱，重叠样本数128，DFT点数2048，采样频率fs title('信号CPSD'); %分析余弦信号的交叉功率谱
fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*t*100)+0.25*randn(size(t));
tau = pi/100; %定义信号相位
y = cos(2*pi*100*(t-tau))+0.25*randn(size(t));
cpsd(x,y,[],[],[],fs) %计算并绘制交叉功率谱密度的幅值
|
【本文地址】
今日新闻 |
推荐新闻 |