[转载]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