程序2 方波信号的富里叶三角级数叠加合成

T0 = 2;A = 2;;%T0为方波的周期;A为方波的幅值;

%周期信号时域描述

tn_i = 1;

for t = -2*T0:0.01:2*T0

if(mod ((t-T0/4),T0) <= T0/2)

y_t(tn_i) = -A; %信号前半周期的表达式

else

y_t(tn_i) = A; %信号后半周期的表达式

end;

t_t(tn_i) = y_t(tn_i);

tn_i = tn_i+1;

end;

%周期信号的频域描述

t = -2*T0:0.01:2*T0;%时域波形的长度

w0 = 2*pi/T0;%基波的频率

figure('color',[1 1 1])

x = 0;%合成的信号值，初始化为0

for i = 1:k

w(i) = (2*i-1)*w0; %第i次谐波的频率

a(i) = (4*A/(pi*((2*i-1)))); %第i次谐波的幅值 注意：a(i)不能为A（i）!

y(i,:) = (-1)^(i+1)*a(i)*cos(w(i)*t); %第i次谐波的值

x = x+y(i,:); %0-i次谐波之和

end;

plot(t,t_t,'b',t,x,'r')

hold on;

xlabel('t');

ylabel('x(t)');

 

 

程序3 信号的功率谱分析

fs = 1000;%采样频率

N = 1001; %点数

t = (0:N-1)/fs;

x1= sin(2*pi*30*t);

x2= 3*sin(2*pi*80*t);

x_n = x1+x2+randn(size(t)); %两个正弦信号频率分别为30Hz、80Hz

[Pxx,w] = periodogram(x_n,[],512,fs); %通过periodogram估计频谱密度

figure('color',[1 1 1])

subplot(2,2,1)%原信号

plot(t,x_n)

xlabel('时间t')

ylabel('被噪声覆盖的X1+X2信号')

subplot(2,2,2)%功率谱密度

plot(w,10*log10(Pxx))

xlabel('频率/Hz')

ylabel('功率谱密度/（dB/Hz）')

subplot(2,2,3)

plot(t,x1)

xlabel('时间t')

ylabel('其中的信号X1')

subplot(2,2,4)

plot(t,x2)

xlabel('时间t')

ylabel('其中的信号X2');

 

程序4 二阶系统的传递函数的单位阶跃及单位脉冲响应的时域波形图

%为简化程序，二阶系统的固有角频率为1,阻尼比分别为0.2、0.7、1

num = [1];

t = 0:0.01:12;

zeta = [0.01 0.2 0.7];

wn = 1;

for k = 1:3

den = [1 2*zeta(k)*wn 1];

sys =tf(num,den);

[y_step(:,k),x] = step(sys,t); %连续系统单位阶跃响应的时域波形函数

[y_impulse(:,k),x] = impulse(sys,t); %单位脉冲响应的时域波形函数

end

subplot(2,1,1)

for k = 1:3

plot (t,y_step(:,k));

hold on

end

xlabel('t');

ylabel('y(t)');

subplot(2,1,2)

for k = 1:3

plot(t,y_impulse(:,k));

hold on

end

xlabel('t');

ylabel('y(t)');

程序5 信号频谱分析和功率谱分析程序

load No1.txt;

fs=2000;%采样频率

N=4096;%抽样数，一般要求取为2^n,如128,256,512,1024,2048,4096,8192等，%应该小于采样数据点数,但至少要包含1个以上波形。

t=No1(:,1)/fs;%从数组No1中取出第一列，再除以采样频率得到采样时间

y=No1(:,2); %从数组No1中取出第二列，得到采样数据

ym=mean(y);%求平均值，找出信号中的直流分量

y=y-ym; %减去平均值，剔除信号中的直流分量

subplot(2,2,1);

plot(1000*t,y) %画时域波形

title('时域波形')

xlabel('时间 (毫秒)')

Yf=fft(y,N);%进行fft变换（离散信号的快速傅里叶变换），求信号的频谱X(f)

f=fs*(0:length(Yf)/2)/length(Yf);%进行对应的频率转换，求出频谱的横坐标

Pyy = Yf.* conj(Yf)/N;%求信号的功率谱

subplot(2,2,2);

plot(f,Pyy(1:(N/2+1))) %画功率谱

[Pymax, fi]=max(Pyy(1:(N/2+1)));%求出最大功率值及其索引号fi

fm=f(fi) %根据索引号求出主频，即功率最大的频率成分的频率

title(['功率谱，主频fm=',num2str(fm),'Hz'])

xlabel('频率 (Hz)')

mag=abs(Yf);%求幅值谱,等于|X(f)|

subplot(2,2,3);

plot(f,mag(1:(N/2+1))) %画幅值谱

[mgmax, mi]=max(mag(1:(N/2+1)));%求出最大功率值及其索引号fi

fm=f(mi) %根据索引号求出主频，即功率最大的频率成分的频率

title(['幅值谱，主频fm=',num2str(fm),'Hz'])

xlabel('频率 (Hz)')

P=angle(Yf); % 求相位谱,等于X(f)的相位角

subplot(2,2,4);

plot(f,P(1:(N/2+1))) %画相位谱

title('相位谱')

xlabel('频率 (Hz)')