西南交通大学信号处理期末作业(3)

plot(t,Sig);

title('原始输入信号'); axis([0 2.1 -7 7]); %% 周期图谱

[Pxx,f]=periodogram(Sig,[],length(t),fs);%周期图法 figure(3); plot(f,Pxx);

title('周期图法求功率谱'); xlabel('f/Hz'); ylabel('功率/db'); %% ARMA谱估计

z=iddata(Sig');%将信号转化为matlab接受的格式

RecordAIC=[];

for p=1:20 %自回归对应PACF，给定滞后长度上限p和q

for q=1:20%移动平均对应ACF m=armax(z(1:length(t)),[p,q]); AIC = aic(m); %armax(p,q)选择对应FPE最小，AIC值最小模型 RecordAIC=[RecordAIC;p q AIC]; end end

for k=1:size(RecordAIC,1) if

RecordAIC(k,3)==min(RecordAIC(:,3)) %选择AIC最小模型

pa_AIC=RecordAIC(k,1); qa_AIC=RecordAIC(k,2); break; end end

mAIC=armax(z(1:length(t)),[pa_AIC,qa_AIC]);

[Pxx2,f2]=freqz(mAIC.c,mAIC.a,fs); P2=(abs(Pxx2).*1).^2; P2tol=10*log10(P2); figure(4);

plot(f2/pi*fs/2,P2tol); title('ARMA法(AIC准则)');xlabel('f/Hz');ylabel('振幅/dB'); plot(RecordAIC(:,3)); ylabel('AIC(p,q)'); %% burg法计算

[Pxx,F] = pburg(Sig,60,length(t),fs);%burg法 figure(6); plot(F,Pxx);

title('Burg法谱估计');

xlabel('f/fs'); %X轴坐标名称 ylabel('功率谱/dB'); %Y轴坐标名称 %%

function y=bandp(x,f1,f3,fsl,fsh,rp,rs,Fs) %带通滤波

%使用注意事项：通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半

%即，f1,f3,fs1,fsh,的值小于 Fs/2 %x:需要带通滤波的序列 % f 1：通带左边界 % f 3：通带右边界

% fs1：衰减截止左边界 % fsh：衰变截止右边界

%rp：边带区衰减DB数设置 %rs：截止区衰减DB数设置 %FS：序列x的采样频率

% f1=300;f3=500;%通带截止频率上下限 % fsl=200;fsh=600;%阻带截止频率上下限 % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值

% Fs=2000;%采样率 %

wp1=2*pi*f1/Fs; wp3=2*pi*f3/Fs; wsl=2*pi*fsl/Fs; wsh=2*pi*fsh/Fs; wp=[wp1 wp3]; ws=[wslwsh]; %

% 设计切比雪夫滤波器；

[n,wn]=cheb1ord(ws/pi,wp/pi,rp,rs); [bz1,az1]=cheby1(n,rp,wp/pi); %查看设计滤波器的曲线 [h,w]=freqz(bz1,az1,256,Fs); h=20*log10(abs(h)); y=filter(bz1,az1,x); end 第5题

%本题目需要提醒一点：给的数据为观测数

据Z而不是X clc; clear;

x1=xlsread('./负荷数据.xls','sheet1'); x1=x1(:,2);

x2=xlsread('./负荷数据.xls','sheet2'); x2=x2(:,2); x=[x1;x2]; N1=length(x1); N=length(x); A=1; B=0; H=1;

w=normrnd(0,1000,1,N);%这里随便取值 v=normrnd(0,1000,1,N); P(1)=16;%随便取值 Z=x;

X(1)=24;%随便取值 R=cov(v); Q=cov(w); for i=2:N

tempx=A*X(i-1);%+B*u(i); TempP=A*P(i-1)*A'+Q;

K(i)=TempP*H'*1/(H*TempP*H'+R); X(i)=X(i-1)+K(i)*(Z(i)-tempx); P(i)=(1-K(i)*H)*TempP; end

t=1:length(Z); figure;

plot(t,Z,'b',t,X(t),'r');title('使用Kalman对电力系统负荷数据进行预测');

xlabel('时间点数');ylabel('电力系统负荷');axis tight;legend('负荷真实值','Kalman预测值'); figure;

subplot(2,1,1);

t=length(x1):length(x);

plot(t,x(t),'b',t,X(t),'r');title('使用Kalman对电力系统负荷数据进行预测');xlabel('时间点数');ylabel('电力系统负荷');axis tight;legend('负荷真实值','Kalman预测值');

set(gca,'XTick',length(x1):2:length(x)); subplot(2,1,2); error=Z-X'; plot(t,error(t));title('预测值与真实值之误差');xlabel('时间点数');

set(gca,'XTick',length(x1):2:length(x));

ylabel('5月5日预测值与真实值误差');axis tight; 第六题：

%% 小波变换 clc; clear; close all;

f=50;%信号频率 oumiga=2*pi*f;

N_sample=2048;%总采样点数 Fs=1000;%采样频率 t=0:1/Fs:1; Tao=0.03;

A=1;%信号幅度 x = 20*exp(-t/Tao)+20*sin(oumiga*t+pi/3)+12*sin(2*oumiga*t+pi/4)+10*sin(3*oumiga*t+pi/6)+6*sin(4*oumiga*t+pi/8)+5*sin(5*oumiga*t+pi/5); % 信号函数表达式 figure; plot(t,x);

title('原始信号');

xlabel('时间t/s','FontSize',14); ylabel('幅值','FontSize',14); %原信号函数

wavename='cmor3-3'; totalscal=256;

Fc=centfrq(wavename); %小波中心频率 c=2*Fc*totalscal; scals=c./(1:totalscal);

f=scal2frq(scals,wavename,1/Fs); % 将尺度转换为频率

coefs=cwt(x,scals,wavename); % 求连续小波系数 figure;

imagesc(t,f,abs(coefs)); colorbar;

xlabel('时间t/s','FontSize',14); ylabel('频率f/Hz','FontSize',14); title('小波时频图','FontSize',16); axis([0 1 0 300]);

%% 短时傅里叶变换

[S,F,T,P]=spectrogram(x,256,250,256,Fs); figure;

surf(T,F,10*log10(P),'edgecolor','none'); axis tight;

view(0,90);

xlabel('时间/s'); ylabel('频率/Hz'); title('短时傅里叶变换结果');

%% Wigner-Ville time-frequency distribution. X=hilbert(x'); [tfr,t,f]=tfrwv(X); figure;

contour(t/Fs,f*Fs,abs(tfr)); xlabel('时间 t/s'); ylabel('频率 f/Hz'); title('Wigner-Ville time-frequency distribution'); axis([0 1 0 300]) %%

第七题： clc; clear; close all; % 参数设置

Fs = 1024; %采样频率 n = 0:1/Fs:2.01;%采样时间 N = length(n); % 采样点

W1=0.001*cos(2*pi*n*10+unifrnd(-pi,pi))+cos(2*pi*50*n+unifrnd(-pi,pi))+0.1*cos(2*pi*n*150+unifrnd(-pi,pi))+0.002*cos(2*pi*n*50+unifrnd(-pi,pi));% 原始信号 x1=awgn(W1,20); %加入噪声 %原信号输出 figure; plot(n,x1);

xlabel('时间(t/秒)','FontSize',10); ylabel('幅值','FontSize',10); axis([0 2.05 -3 3]); title('原始信号'); %% 小波变换

wavename='cmor3-3'; totalscal=256;

Fc=centfrq(wavename); %小波中心频率 c=2*Fc*totalscal; scals=c./(1:totalscal);

f=scal2frq(scals,wavename,1/Fs); % 将尺度转换为频率

coefs=cwt(x1,scals,wavename); % 求连续小波系数 figure;

imagesc(n,f,abs(coefs)); colorbar;

xlabel('时间t/s','FontSize',14); ylabel('频率f/Hz','FontSize',14); title('小波时频图','FontSize',16); %% 短时傅里叶变换

[S,F,T,P]=spectrogram(x1,256,250,256,Fs); figure;

surf(T,F,10*log10(P),'edgecolor','none'); axis tight;

view(0,90);

xlabel('时间/s'); ylabel('频率/Hz'); title('短时傅里叶变换结果'); %% 维格纳威利分布 X=hilbert(x1'); [tfr,t,f]=tfrwv(X); figure;

contour(t/Fs,f*Fs,abs(tfr)); xlabel('时间 t/s'); ylabel('频率 f/Hz');

title('Wigner-Villetime-frequency distribution');

%% 周期图谱估计

[Pxx,f]=periodogram(x1,[],length(x1),Fs);%周期图法 figure; plot(f,Pxx);

title('周期图法求功率谱'); xlabel('f/Hz'); ylabel('功率/db'); set(gca,'XTick',0:50:600); %% MUSIC方法谱估计 nfft=1024; figure;

pmusic(x1,[7,1.1],nfft,Fs,32,16); grid on;

xlabel('频率(f/Hz)','FontSize',10); ylabel('功率(dB)','FontSize',10); title('MUSIC方法');

%% burg法谱估计

[Pxx,F] = pburg(x1,60,length(x1),Fs);%burg法 figure;

plot(F,Pxx);

title('Burg法谱估计');

xlabel('f/fs'); %X轴坐标名称

ylabel('功率谱/dB'); %Y轴坐标名称 set(gca,'XTick',0:50:600);

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说教育文库西南交通大学信号处理期末作业(3)在线全文阅读。