%对X(t)进行频谱分析
Fx=fft(x,N); %对x(t)进行fft变换，在0~16000区间内得到2N-1个频率值 magn=abs(Fx); %求x(t)幅值 xangle=angle(Fx); %求X(t)相位
labelang=(0:length(x)-1)*16000/length(x); %在0~16000区间内求横坐标刻度
figure(2); plot(labelang,magn*10); %在0~16000区间内做频谱和相位图 axis([0 16000 -0.5 600]); xlabel('频率/Hz');ylabel('幅值');title('X(t)频谱图');grid on;
%求X(t)的自相关函数
[c,lags]=xcorr(x,'coeff'); %求自相关函数
figure(3); subplot(2,1,1); plot(lags/fs,c); %绘制自相关函数 axis tight; xlabel('T');ylabel('Rx(T)');title('X(t)的自相关函数');grid on; %求X(t)的功率谱密度 long=length(c);
Sx=fft(c,long); labelx=(0:long-1)*2*pi;
plot_magn=10*log10(abs(Sx));
subplot(2,1,2); plot(labelx,plot_magn); %绘制功率谱密度
axis tight;xlabel('w');ylabel('Sx(w)');title('X(t)的功率谱密度');grid on; %窄带系统检测
z1=2.*cos(2*pi*4000*t); z2=-2.*sin(2*pi*4000*t); Ac=z1.*x; ?分量 As=z2.*x; %As分量
y=Ac.*cos(2*pi*4000*t)-As.*sin(2*pi*4000*t); %滤波器设计
f_p=1000;f_s=1600;R_p=1;R_s=35; %设定滤波器参数; 通、阻带截止频率,通、阻带衰减
Ws=2*f_s/fs;Wp=2*f_p/fs; %频率归一化
[n,Wn]=buttord(Wp,Ws,R_p,R_s); %采用巴特沃思滤波器 [b,a]=butter(n,Wn); %求得滤波器传输函数的多项式系数 figure(4);
[H,W]=freqz(b,a); %求得滤波器传输函数的幅频特性
subplot(2,1,1); plot(W*fs/(2*pi),abs(H)); %在0～2pi区间内作幅度谱 title('低通滤波器幅度谱'); grid on;
subplot(2,1,2); plot(W*fs/(2*pi),angle(H)); %在0～2pi区间内作相位谱 title('低通滤波器相位谱'); grid on; %求Ac(t)滤波后的统计特性
mc=filter(b,a,Ac); %上支路通过滤波器 Ac(t) disp('Ac(t)的均值');Eh=mean(mc) %求均值
disp('Ac(t)的均方值是');E2h=mc*mc'/N %求均方值 disp('Ac(t)的方差');Dh=var(mc) %求方差 %画Ac(t)的 波形
21
figure(6); subplot(2,1,1); n=0:N-1; plot(n,mc);
axis([0 300 -1 1]);xlabel('采样点');ylabel('幅值');title('Ac(t)的时域波形');grid on;
%画Ac(t)的 图
yc=fft(mc,length(mc)); %作傅里叶变换
longc=length(yc); %求傅里叶变换后的序列长度 labelx=(0:longc-1)*16000/longc; magnl=abs(yc); %取幅值
subplot(2,1,2); plot(labelx,magnl); %绘制幅频曲线
axis tight; xlabel('频率(Hz)'); ylabel('幅值'); title('Ac(t)频谱图'); grid on;
%求Ac(t)的 函数
[c1,lags1]=xcorr(mc,'coeff'); %求自相关函数
figure(7); subplot(2,1,1); plot(lags1/fs,c1); %绘制自相关函数曲线 xlabel('T');ylabel('Rx(T)');axis tight; title('Ac(t)的自相关函数'); grid on;
%求Ac(t)的 功率谱
Sac=fft(c1,length(c1)); % 求频率响应 magnc=abs(Sac); % 求幅频响应 long=length(Sac); % 频谱序列长度 labelc=(0:long-1)*16000/long;
subplot(2,1,2); plot(labelc,10*log10(magnc)); % 绘制双边功率谱
xlabel('频率(Hz)');ylabel('功率谱(dbW)');axis tight;title('Ac(t)的双边功率谱');grid on;
%求得As(t)的统计特性
ms=filter(b,a,As); %对下支路信号进行滤波得As(t) disp('As(t)的均值'); Eh=mean(ms) % 求均值
disp('As(t)的均方值是'); E2h=ms*ms'/N % 求均方值 disp('As(t)的方差'); Dh=var(ms) % 求方差 %作As(t)的时域波形
figure(8);subplot(2,1,1); n=0:N-1;plot(n,ms); % 绘制时域波形
axis([0 300 -0.5 2]); xlabel('采样点');ylabel('幅值');title('As(t)的时域波形');grid on;
%对As(t)进行FFT变换并做频谱图
ys=fft(ms,length(ms)); % 进行傅里叶变换
longs=length(ys); % 求傅里叶变换后的序列长度 labelx=(0:longs-1)*16000/longs; magn2=abs(ys); % 取幅值
subplot(2,1,2); plot(labelx,magn2); % 绘制幅频曲线
axis tight;xlabel('频率(Hz)');ylabel('幅值');title('As(t)的频谱图');grid on; %求As(t)的自相关函数
[c2,lags2]=xcorr(ms,'coeff'); %求出As(t)的自相关序列
figure(9);subplot(2,1,1);plot(lags2/fs,c2); %绘制自相关函数曲线
22
xlabel('T');ylabel('Rx(T)');axis tight;title('As(t)的自相关函数');grid on; %求As(t)的双边功率谱
Sas=fft(c2,length(c2)); %对As(t)的自相关函数进行傅里叶变换 magnc=abs(Sac); %求As(t)的双边功率谱幅值 long=length(Sas); %求傅里叶变换后的序列长度 labels=(0:long-1)*16000/long;
subplot(2,1,2); plot(labelc,10*log10(magnc)); %画As(t)的自相关函数频谱
xlabel('频率(Hz)');ylabel('功率谱(dbW)');axis tight;title('As(t)的双边功率谱'); % 求y(t)的统计特性
disp('输出信号Y(t)的均值');Eh=mean(y) % 求均值 disp('输出信号Y(t)的均方值');E2h=y*y'/N %求均方值 disp('输出信号Y(t)的方差');Dh=var(y) % 求方差 %作输出信号Y(t)的时域波形
figure(10); subplot(2,1,1);n=0:N-1;plot(n,y);
axis([0 150 -2 2]);xlabel('采样点');ylabel('幅值');title('Y(t)的时域波形');grid on;
%进行FFT变换并做频谱图
yy=fft(y,length(y)); %进行傅里叶变换
longy=length(yy); %求傅里叶变换后的序列长度 labelx=(0:longy-1)*16000/longy; magn3=abs(yy); %取幅值
subplot(2,1,2); plot(labelx,magn3); %绘制幅频曲线
axis tight;xlabel('频率(Hz)');ylabel('幅值');title('Y(t)的频谱图');grid on; %求输出信号Y(t)的自相关函数
[c3,lags3]=xcorr(y,'coeff'); %求自相关函数
figure(11); subplot(2,1,1); plot(lags3/fs,c3); %绘制自相关函数曲线 xlabel('T');ylabel('Rx(T)');axis tight;title('Y(t)的自相关函数');grid on; %求输出信号Y(t)的双边功率谱
Sy=fft(c3,length(c3)); %进行傅里叶变换 magny=abs(Sy); %取幅值 long=length(Sy);
labely=(0:long-1)*16000/long;
subplot(2,1,2); plot(labely,10*log10(magny)); %****画Y(t)的功率谱密度 xlabel('频率(Hz)');ylabel('功率谱(dbW)');axis tight;title('Y(t)的双边功率谱');grid on;
23
结果
24
25