[转载]MATLAB对信号的时域与频域分析
abel('k');ylabel('spectrum');
title('DFS:L=12,N=60');
%频谱分析的若干问题
%不同时窗正弦信号频谱的分析
%1Hz的正弦波仅取10个周期及50个周期,都做4096点FFT,比较两者的频谱
t=0:0.1:10; %取10或50(将10变成50即可)个周期的进行变换
y1=sin(2*pi*t);
Y1=fft(y1,4096); %4096点的快速傅里叶变换
Y=fftshift(Y1);
c=[0:2047]/409;
plot(c,abs(Y(2049:4096)))
axis([0 2 -5 60]);
title('sin(t)周期 T=10')
grid
%加窗函数后的频谱特征分析
%不加窗函数是正弦谱线与矩形窗频谱的卷积;加窗函数是正弦谱线于汉明窗频谱的卷积
%后者的旁瓣小得多matlab离散信号的时域和频域分析,但是主谱线宽一些
t=0.1:0.1:20; %取20个周期
y=sin(2*pi*t); %定义1Hz的正弦波
w=hamming(200); %定义200长的汉明窗
y=y.*w';
Y1=fft(y,4096); %4096点的快速傅里叶变换
Y=fftshift(Y1);
c=[0:2047]./409.6;
plot(c,abs(Y(2049:4096)))
axis([0 3 -5 60]);title('sin(t) 加窗频谱');grid
%CZT线性调频Z变换
%应用CZT变换不在频率轴上均匀采样,对关心的部分加大采样密度,可以提高关心部分的
%频率matlab离散信号的时域和频域分析,可将两个频率相距很近的信号区分开来
%主要指令:CZT(Chirp-Z变换)
N=256;
f1=13.2;f2=13.41;f3=15;fs=60;
n=0:N-1;
t=2*pi*n/fs;
e=fs/N;
n1=0:e:(fs/2)-e;
x=sin(f1*t)+sin(f2*t)+sin(f3*t); %三个不同频率正弦信号相加
Y1=abs(fft(x)); %一般傅里叶频谱分析
subplot(121)
plot(n1,(Y1(1:N/2)));
title('fft');grid;
M=60; %CZT变换的长度
f0=12.6;
q=0.05;
A=exp(j*2*pi*f0/fs); %CZT变换的起点
W=exp(-j*2*pi*q/fs); %CZT变换的倾斜率
Y3=czt(x,M,W,A); %CZT变换做频谱分析
n2=f0:q:f0+(M-1)*q;
subplot(122)
plot(n2,abs(Y3));
title('czt');grid
%常见信号的频谱分析
%基础波形 fft(快速傅里叶变换),fftshift(移动傅里叶变换的系数)
n=2048;t=1:2048;
y=zeros(1,n);
y(1,[1020:1035])=1; %定义宽度为15的方波
y1=[zeros(1,1000),ones(1,60),zeros(1,988)];%定义宽度为60的方波
subplot(121)
plot(t,y,'b',t,y1,'r')
axis([850 1200 -0.2 1.3]);
title('方形图');grid
subplot(122)
Y2=fft(y,2048);
Y2a=fftshift(Y2);
Y3=fft(y1,2048);
Y3a=fftshift(Y3);
w=-1024:1023;
semilogy(w,abs(Y2a),'b',w,abs(Y3a),'r')
axis([-500 500 -5 80]);
title('方波频谱');
axis([-250 250 1 100]);grid;
%三角波
n=2048;t=1:2048;
y=[zeros(1,800),[1:500],zeros(1,748)]; %定义三角波
subplot(121)
plot(t,y)
axis([0 2100 -30 560]);
title('三角波');grid
subplot(122)
Y1=fft(y,2048);
Y=fftshift(Y1);
c=-1024:1023;
yw=log(abs(Y)+eps);
plot(c,yw,'r')
axis([-104 104 7 12]);
title('三角波频谱');
grid on;
%冲击函数
n=400;
delta=4*pi/n;
t=-2*pi:delta:2*pi;
y=sinc(t); %定义冲击函数
subplot(121)
plot(t,y)
axis([-7 7 -0.4 1.3]);
title('sinc(t)');grid
subplot(122)
Y1=fft(y,2048);
Y=fftshift(Y1);
c=-1024:1023;
plot(c,abs(Y))
axis([-500 500 -5 40]);
title('sinc(t)的频谱');grid
%调制信号
%双边带调幅
fs=1000;
t=0:1/fs:.4;
fc=250;
y=0.4*cos(pi*20*t);
x=modulate(y,fc,fs,'amdsb-tc'); %双边带调幅
subplot(121)
plot(t,x)
legend('调幅双边带波形');grid
subplot(122)
p=fft(x,1024);
p1=fftshift(p);
w=0:511;
p2=abs(p1);
plot(w,p2(1:512));
grid;
legend('调幅双边带频谱')
axis([100 400 -10 100])
%用直方图表示调频信号谱
m=5; %用直方图表示调制指数为5的调频信号谱分析
n=1:10;
y1(n)=besseli(n,(i*m),1);
z1=real(y1./(i.^(n)));
z2=fliplr(z1);
y0=besseli(0,(i*m),1);
z0=real(y0);
z=[z2 z0 z1];
x=-10:10;
bar(x,abs(z))
axis([-10 10 -.1 .8]);grid
title('m=5')
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%四种数字调制方法(ASK FSK PSK OQPSK)的频谱图
%产生码元宽度为64的随机序列
%m文件
n=1:8192;
m=1:128;x(n)=randint(1,8192,2);x=[x(n)]';
y(n)=zeros(1,8192);z(m)=zeros(1,128);
for n=1:8192
for m=1:128
if n==64*m-63
z(m)=x(n);
if m==ceil(n/64)
y([(64*m-63):(64*m)]')=z(m);
end
end
end
end
n=1:8192;rm2=y(n);
%ASK调制
n=[1:(2^13)];
x1=cos(n.*1e9*2*pi/4e9);
run('rm2');
X2=rm2;x2=X2';
x=x1.*X2;
b=.42+.5*cos(2*pi*(n-(2^12))/(2^13))+0.08*...
cos(4*pi*(n-(2^12))/(2^13));
X=b.*x;x3=[ones(1,64) zeros(1,8128)];
y1=X(1:(2^13));y4=x1.*x3;
Y1=fft(y1,(2^13));magY1=abs(Y1(1:1:(2^12)+1))/(200);
Y4=fft(y4,(2^13));magY4=abs(Y4(1:1:(2^12)+1))/(40);
k1=0:(2^12);w1=(2*pi/(2^13))*k1;
u=(2*w1/pi)*1e9;
figure(1)
subplot(211)
plot(u,magY1,'b',u,magY4,'r');grid
title('ASKr');axis([4e8,1.6e9,0,1.1])
X2=b.*X2;
y2=X2(1:(2^13));
Y2=fft(y2,(2^13));magY2=abs(Y2(1:1:(2^12)+1))/(200);
k1=0:(2^12);w1=(2*pi/(2^13))*k1;
u=(2*w1/pi)*1e9;
Y3=fft(x3,(2^13));magY3=abs(Y3(1:1:(2^12)+1))/(40);
subplot(212)
semilogy(u,magY2,'b',u,magY3,'r');grid
title('ASKr-modulation');
axis([0 1.2e9 3e-4 3]);
figure(2)
subplot(211);plot(n,x2);axis([0 720 -0.2 1.2])
subplot(212);plot(n,x);axis([0 720 -1.2 1.2])
%FSK调制频谱分析
n=[1:(2^13)];d=.13;
run('rm2');
x2=((2*rm2)-1);
x5=cos(.5*n.*((1e8*2*pi/4e8)-d));
x4=cos(.5*n.*((1e8*2*pi/4e8)+d));
x=cos(.5*n.*((1e8*2*pi/4e8)+d*x2));
b=.42+.5*cos(2*pi*(n-(2^12))/(2^13))+0.08*...
cos(4*pi*(n-(2^12))/(2^13));
X=b.*x;x3=[ones(1,64) zeros(1,8128)];
y1=X(1:(2^13));y4=x4.*x3;y5=x5.*x3;
Y1=fft(y1,(2^13));magY1=abs(Y1(1:1:(2^12)+1))/(200*2);
Y4=fft(y4,(2^13));magY4=abs(Y4(1:1:(2^12)+1))/(70);
Y5=fft(y5,(2^13));magY5=abs(Y5(1:1:(2^12)+1))/(70);
k1=0:(2^12);w1=(2*pi/(2^13))*k1;
u=(2*w1/pi)*1e9;
figure(1)
subplot(211)
plot(u,magY1,'b',u,magY4,'r',u,magY5,'g');grid
title('FSKr');axis([1e8 1e9 0 .65]);
X2=b.*x2;
y2=X2(1:(2^13));
Y2=fft(y2,(2^13));magY2=abs(Y2(1:1:(2^12)+1))/(50*6);
k1=0:(2^12);w1=(2*pi/(2^13))*k1;
u=(2*w1/pi)*1e9;
Y3=fft(x3,(2^13));magY3=abs(Y3(1:1:(2^12)+1))/(30);
subplot(212)
plot(u,magY2,'b',u,magY3,'r');grid
title('FSKr-modulation');axis([0 1.2e9 0 2.3])
figure(2)
subplot(211);plot(n,x2);axis([0 320 -1.2 1.2])
subplot(212);plot(n,x);axis([0 320 -1.2 1.2])
%PSK调制
n=[1:(2^13)];
x1=cos(n.*1e9*2*pi/4e9);
run('rm2');
x2=(2*rm2)-1;
x=x1.*x2;
b=.42+.5*cos(2*pi*(n-(2^12))/(2^13))+0.08*...
cos(4*pi*(n-(2^12))/(2^13));
X=b.*x;x3=[ones(1,64) zeros(1,8128)];
y1=X(1:(2^13));y4=x1.*x3;
Y1=fft(y1,(2^13));magY1=abs(Y1(1:1:(2^12)+1))/(300);
Y4=fft(y4,(2^13));magY4=abs(Y4(1:1:(2^12)+1))/(30);
k1=0:(2^12);w1=(2*pi/(2^13))*k1;
u=(2*w1/pi)*1e9;
figure(1)
subplot(211)
plot(u,magY1,'b',u,magY4,'r');grid
title('PSKr');axis([4e8 1.6e9 0 1.1]);
X2=b.*x2;
y2=X2(1:(2^13));
Y2=fft(y2,(2^13));magY2=abs(Y2(1:1:(2^12)+1))/(300);
k1=0:(2^12);w1=(2*pi/(2^13))*k1;
u=(2*w1/pi)*1e9;
Y3=fft(x3,(2^13));magY3=abs(Y3(1:1:(2^12)+1))/(30);
subplot(212)
plot(u,magY2,'b',u,magY3,'r');grid
title('PSKr-modulation');axis([0 1.2e9 0 2.3])
figure(2)
subplot(211);plot(n,x2);title('PSKr');
axis([0 1200 -1.2 1.2])
subplot(212);plot(n,x);axis([100 275 -1.2 1.2])
%谱估计 tfe 从输入输出中估计传递函数
%估计传输函数
h=fir1(30,.2,boxcar(31)); %使用矩形窗
x=randn(16384,1); %输入信号
y=filter(h,1,x); %输出信号
z=tfe(x,y,1024,[],[],512);
n=[1:length(z)];
p=log10(abs(z)+eps);
plot(n,p);grid
%互功率谱密度 csd已知二信号求互功率谱密度
h=fir1(30,.2,hamming(31)); %使用汉明窗
h1=ones(1,10)/sqrt(10);
r=randn(16384,1);
x=filter(h1,1,r);
y=filter(h,1,x);
z=csd(x,y,1024,10000,triang(500),0,[]); %使用三角窗
n=[1:length(z)];
p=log10(abs(z)+eps);
plot(n,p);grid
%功率谱密度 psd求信号的功率谱密度
h=fir1(30,.2,triang(31)); %使用三角窗
x=randn(16384,1);
y=filter(h,1,x);
z=psd(y,1024,10000,kaiser(512,5),0,.95); %使用凯塞窗
n=[1:length(z)];
p=log10(abs(z)+eps);
plot(n,p);grid
马云创业时