我这个matlab程序能不能运行，有没有拼写错误，，

%谐波叠加法模咐厅销拟风速时程（采用Kaimal谱）

clc

clear

%风速时程参数设定

m=10%模拟风速时程个数

N=2^8 %频率采样点数

dt=0.5 %时间间隔

omegaup=2*pi%上限频率

%风速谱参数设定

L=1000 %跨度

z=50%模拟点离地面高度(单位伏粗 米）

z0=0.03 %地面粗糙度

Uz=40 %50米处平均风速，摩阻速度（单位 米/秒）

delta=100 %模拟点间距（单位 米）

lamda=10%空间相关函数系数

K=0.4 %Kaman常数0.4

M=2*N

v=zeros(m,M*m)

t=0.5*(0:1:(M*m-1))

domega=omegaup/N

D=zeros(m,m,N)

Coh=zeros(m,m,N)

S=zeros(m,m,N)

H=zeros(m,m,N)

U=K*Uz/log(z/z0) %U为摩阻速度(通过50米的平均风速和z求得参数u*摩阻）

%形成目标谱

omega1=omegaup/N:omegaup/N:omegaup

Sw1=200*U^2*z/Uz./(1+50*omega1*z/(2*pi*Uz)).^(5/3)%Kaman谱的表达式S(w)功率谱密度

for j=1:m % 对m条风速循环

%rand('state'衡游,0) %模拟随机数序列

thet=2*pi*rand(j,N) %生成随机相位

for l=1:N %生成频率序列

omega1(l)=(l-1)*domega+j/m*domega

end

Sw=200*U^2*z/Uz./(1+50.*omega1.*z./(2*pi*Uz)).^(5/3)

%计算谱数据库矩阵，功率谱计算，Kaimal谱

for ji=1:m

for l=1:m

for k=1:N

Coh(ji,l,k)=(exp(-lamda*omega1(k)*delta/(2*pi*Uz)))^(abs(ji-l))

S(ji,l,k)=Sw(k)*Coh(ji,l,k)

end

end

end

%Cholesky 矩阵分解

for i=1:1:N

H(:,:,i)=chol(S(:,:,i))

H(:,:,i)=H(:,:,i)'

end

% 填充谱数据矩阵D

D(:,j,:)=H(:,j,:)

i=sqrt(-1)

B1=sqrt(2*domega).*D(j,:,:)

B2=zeros(j,N)

for ii=1:j

for jj=1:N

B2(ii,jj)=B1(1,ii,jj)

end

end

B2=B2.*exp(i.*thet)

G=zeros(j,M)

for jj=1:j

G(jj,1:M)=fft(B2(jj,:),M)

for jjj=2:m

G(jj,((jjj-1)*M+1):(jjj*M))=G(jj,1:M)

end

end

%谐波叠加生成模拟点的风速时程

for p=1:M*m

for k=1:j

v(j,p)=v(j,p)+real(G(k,p)*exp(i*k/m*domega*(p-1)*dt))

end

end

end

%显示风速时程

figure(1) %第1点风速时程

plot(t,v(1,:))

figure(2) %第5点风速时程

plot(t,v(5,:))

[power1,frq1]=spectrum(v(1,:),M*m,2,boxcar(512),0,'mean')%计算第1点的模拟自功率谱

[power5,frq5]=spectrum(v(5,:),M*m,2,boxcar(512),0,'mean')%计算第5点的模拟自功率谱

%计算第一。第五点的模拟互功率谱

%****不太懂你下边的程序，你自己修改一下吧，上边的已经debug过了****%

for i=1:N

s15=S(1,5,i)

Sw15=Sw15(s15)

end

%功率谱检验

figure(3) %第1点模拟自功率谱和计算自功率谱比较（对数坐标）

loglog(freq1,power1,'r',omega1,Sw1,'b')

figure(4) %第5点模拟自功率谱和计算自功率谱比较（对数坐标）

loglog(freq5,power5,'r',omega1,Sw1,'b')

figure(5) %第1,5点模拟自功率谱和计算自功率谱比较（对数坐标）

plot(ferq15,abs(power15),'r',omega,Sw15,'b')

先说分解：对三角波求一次导是迅斗前方波销悔信号，求二次倒是冲激信号……图会画吧？

再说合成，根据网络提供的程序，懒得自己写了：

A0=5

f0=20

N=2048

fs=1024

P0=pi/2%因为三角波的傅里叶展开式是余弦函数，欲用MySin叠加生成，则需将初相位设为90°

b=f0%将基频赋给b

a=1:2:1024

PanSin=[]%正弦信号叠加数组

for n=1:30 %n为叠加个数，设最大30

f0=a(n)*b%各个谐波频率

[t,xt]=MySin(A0,f0,P0,N,fs)%调用[t,xt]=MySin(A0,f0,P0)函数

PanSin(1:length(t),n)=1/(a(n)^2)*xt %进行叠加

end %30次叠加完成后，PanSin数组中是30列长度为length(t)的数据

b=PanSin'%求转置

y2=A0/2+sum(b(1:2,:))%取前2行叠加

y5=A0/2+sum(b(1:5,:))%取前5行叠加

y10=A0/2+sum(b(1:10,:))%取前10行叠加

y20=A0/2+sum(b(1:20,:))%取前20行叠加

y30=A0/2+sum(PanSin')%取前30行叠加

plot(t(1:103),y2(1:103),'-r')%绘图，只显示2个周期

hold on

plot(t(1:103),y5(1:103),'-g')

plot(t(1:103),y10(1:103),'-y')

plot(t(1:103),y20(1:103),'-m')

plot(t(1:103),y30(1:103),'-b')

title(['用正弦信号构造三角波信号'],'fontsize',14)

xlabel('t(ms)','fontsize',14)%设置x轴标题为t(ms)

ylabel('x(t)','fontsize',14)%设置y轴标题为x(t)

axis tight

legend('亩清2个正弦信号','5个正弦信号','10个正弦信号','20个正弦信号','30个正弦信号')

这个东西网上很多哦。自己先琢磨一下，实在不懂再问我吧。

function rectexpd(T1,T0,m)

%矩形信号串信号分解与合成

%T1：矩形信号区间为（-T1/2，T1/2）

%T0：矩形信号串周期

%m:傅里叶级数档罩巧展开项次数

t1=-T1/2:0.01:T1/2

t2=T1/2:0.01:(T0-T1/2)

t=[(t1-T0)'(t2-T0)'t1't2'(t1+T0)']

n1=length(t1)n2=length(t2)%根据周期矩形信号函数周期,计算点数

f=[ones(n1,1)zeros(n2,1)ones(n1,1)zeros(n2,1)ones(n1,1)] %构造周期矩形信号串

y=zeros(m+1,length(t))y(m+1,:)=f'

figure(1)

plot(t,y(m+1,:))%绘制周期矩形信号串

axis([-(T0+T1/2)-0.5,(T0+T1/2)+0.5,0,1.2])

set(gca,'XTick',[-T0,-T1/2,T1/2,T0])

set(gca,'XTickLabel',{'-T0','-T1/2','T1/2','T0'})

title('矩形信号串')grid on

a=T1/闷谨T0

pause%绘制离散幅度谱

freq=[-20:1:20]

mag=abs(a*sinc(a*freq))

stem(freq,mag)

x=a*ones(size(t)) %T1f0

for k=1:m %循环显示谐波叠加图形

pause

x=x+2*a*sinc(a*k)*cos(2*pi*t*k/T0)

y(k,:)=x

plot(t,y(m+1,:))hold on

plot(t,y(k,:))hold offgrid on

axis([-(T0+T1/2)-0.5,[T0+T1/2]+0.5,-0.5,1.5])

title(strcat(num2str(k),'次谐波叠加'))

xlabel('t')end %同一级对齐

pause

plot(t,y(1:m+1,:))grid on

axis([-T0/2,T0/2,-0.5,1.5])

title('行键各次谐波叠加')xlabel('t')

---