日志
功率谱转时域有关的两个程序
|
%%本程序是根据白噪声的功率谱来求时域信号
%%源程序来自于:王济,胡晓.MATLAB在振动信号处理中的应用.中国水利水电出版社,2005,150:152
clear
clc
sf=400;%采样频率
fi=0.5;%最小截止频率
fa=90;%最大截止频率
tl=40;%时间长度
am=1.1;%最大幅值
n=round(sf*tl)+1;%计算白噪声信号数据长度
t=0:1/sf:(n-1)/sf;%建立信号离散时间向量
nfft=2^nextpow2(n);%确定FFT长度
ni=round(fi*nfft/sf+1);%确定最小截止频率对应数组元素的下标
na=round(fa*nfft/sf+1);%确定最大截止频率对应数组元素的下标
r=zeros(1,nfft/2);%建立一个长度为nfft/2元素全为0的向量
r(ni:na)=1;%将r的正频率带通内的元素赋为1,其余为0 *****************(1)
a=2*pi*rand(1,nfft/2);%产生随机向量
b=r.*exp(i*a)*nfft/2;%将幅值谱和相位谱转换在实部和虚部
a=[b,b(nfft/2:-1:1)];%将正负圆频率向量组合成一个向量
x=real(ifft(a,nfft));%IFFT变换,并取实部为白噪声信号
y=am*x(1:n)/max(abs(x(1:n)));%取指定长度的数据并使数据的最大值为指定的幅值 ********(3)
nft=1024;%定义自功率谱的FF长度
y1=psd(y,nft,sf,boxcar(nft),nft/2);%计算白噪声的自功率谱**********(2)
f=0:sf/nft:(nft/2-1)*sf/nft;%确定频率向量
subplot(2,1,1)
plot(t,y,'k-')
xlabel('频率 (Hz)')
ylabel('幅值')
grid on
subplot(2,1,2)
plot(f(1:nft/4),y1(1:nft/4),'k-')
xlabel('频率 (Hz)')
ylabel('幅值')
grid on
fs=300;%设定采样频率
N=1000;
n=0:N-1;
y=sin(-i*0.0013*x)
figure(1)
subplot(312)
plot(t,y)%
fmin=0;%最小截止频率
fmax=1;%最大截止频率
ny=length(y);
nfft=2^nextpow2(ny);
%计算频率间隔（Hz/s）
df=fs/nfft;
%计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
%计算圆频率间隔（rad/s）
dw=2*pi*df;
%建立正的离散圆频率向量
w1=0:dw:2*pi*0.5*fs;
%建立负的离散圆频率向量
w2=-2*pi*(0.5*fs-df):dw:-dw;
%将正负圆频率向量组合成一个向量
w=[w1,w2];
Y=fft(y,nfft)
%进行积分的频域变换
a=zeros(1,nfft);
for i=2:nfft-1
a(i) =Y(i)./w(i);
end
a1=imag(a);
a2=real(a);
Y1=a1-a2*i;
a=zeros(1,nfft);
%消除指定正频带外的频率成分
a(ni:na)=Y1(ni:na);
%消除指定负频带外的频率成分
a(nfft-na+1:nfft-ni+1)=Y1(nfft-na+1:nfft-ni+1);
Y2=ifft(a,nfft);
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
Y3=real(Y2(1:N))
figure(1)
subplot(313)
plot(t,Y3);
