%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%螺旋桨轴频和唱音谐波关系的线谱程序，2015年6月6日。机械噪声的线谱也可以用这种正弦信号叠加产生
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%唱音与螺旋桨转速有关系，由公式可以计算出转速，但是公式编写程序后与论文上给的结果差100倍
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DT=1000;
%v=5;
%kexi=(1/(6.485936*10^(-3)+2.479736*10^(-5)*DT+5.046719*10^(-11)*DT^2))^(2/3);
%s=v/kexi*DT*10^(-4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%产生螺旋桨唱音的线频噪声的相应频率,设定幅值和初始相位
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=4; %螺旋桨叶片数，4叶桨
for m=1:9 % m谐频次数
s=[3,5,7]; % s是轴频
for i=1:3 %产生3种航速低速、中速和高速，用来生成fm相应螺旋桨唱音的频率
fi(m,i)=s(i).*m*n; %计算fm=mns
end
end
A=[1.30,1.40,1.60
1.25,1.35,1.55
1.20,1.30,1.50
1.15,1.25,1.40
1.10,1.20,0
1.05,1.17,0
1.00,0,0
0.95,0,0
0.90,0,0]; %幅值
a=-pi;
b=pi;
fai=a+(b-a)*rand(9,3); %取fai的［-pi,pi］的随机值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算sin(2pifit+fai)信号的和，用正弦信号做线谱仿真
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fs=1000;
N=500;
n=linspace(1,1000,N);
y1=1.00*sin(2*pi*3*n/fs);%低航速轴频正弦信号可以求出轴频线谱
y2=1.20*sin(2*pi*5*n/fs);%中航速轴频正弦信号可以求出轴频线谱
y3=1.40*sin(2*pi*7*n/fs);%高航速轴频正弦信号可以求出轴频线谱
Xline=zeros(); %用来存储唱音频率的正弦信号
Linex=zeros();%用来存储每一航速情况下，9个唱音频率，时间T采样9次的正弦信号值
Line=zeros();%用来存储仿真线谱的正弦信号和，3行分别代表低速、中速和高速，9列代表时间T采样9次
Sumline=zeros();%用来计算存储一组航速下的正弦信号的和
for i=1:3 % 三种航速，列下标
for j=1:9
Xline(j,1:N)=A(j,i)*sin(2*pi*fi(j,i)*n/fs+fai(j,i));%分别求9个唱音频率的正弦信号的值
Sumline(j+1,1:N)=Xline(j,1:N);
for k=1:N
Sumline(j+1,k)=Sumline(j,k)+ Sumline(j+1,k);%求每个航速下的线谱的正弦信号和
end
end
Linex((i-1)*9+1:i*9,1:N)= Xline(1:9,1:N); %存储每一航速下的9种唱音频率正弦信号值
Line(i,1:N)=Sumline(10,1:N); %第一行代表低速，第二行代表中速，第三种低表高速，
% 假设低速时螺旋桨不产生空化，高速时中螺旋桨产生空化，中速时介于空化与非空化之间
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%sin(2pifit+fai)图形正弦时域显示
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
subplot(3,1,1)
plot(n/fs,Line(1,:))
xlabel('t');
ylabel('幅值');
title('低速时域波形');
hold on
subplot(3,1,2)
plot(n/fs,Line(2,:))
xlabel('t');
ylabel('幅值');
title('中速时域波形');
hold on
subplot(3,1,3)
plot(n/fs,Line(3,:))
xlabel('t');
ylabel('幅值');
title('高速时域波形');
hold on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求螺旋桨唱音的频谱，显示3种航速下的轴频和唱音线谱的频谱图形
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% x=sin(2*pi*f0*t);
% % y=sin(2*pi*10*n/fs);
% ft1=fft(y,N);
ft1=fft( y1+ Line(1,1:N),N);%进行fft变换
a1=abs(ft1);%求幅值
%a1=a1*2/N; %幅值单位转换为dB
f=(0:length(ft1)-1)*fs/length(ft1);%进行对应的频率转换
figure(2);
subplot(311);
plot(f,a1);%做频谱图
xlabel('频率/Hz');
ylabel('幅值dB/');
title('低速线谱谱图');
grid;
ft2=fft(y2+Line(2,1:N),N);%进行fft变换
f=(0:length(ft2)-1)*fs/length(ft2);%进行对应的频率转换
a2=abs(ft2);%求幅值
%a2=a2*2/N; %幅值单位转换为dB
subplot(312);
plot(f,a2);%做频谱图
xlabel('频率/Hz');
ylabel('幅值dB/');
title('中速线谱谱图');
grid;
ft3=fft(y3+Line(3,1:N),N);%进行fft变换
f=(0:length(ft3)-1)*fs/length(ft3);%进行对应的频率转换
a3=abs(ft3);%求幅值
%a3=a3*2/N; %幅值单位转换为dB
subplot(313);
plot(f,a3);%做频谱图
% axis([0,100,0,80]);
xlabel('频率/Hz');
ylabel('幅值dB/');
title('高速线谱谱图');
grid;
评论0