- %生成白噪声信号
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- clear
- clc
- close all hidden
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % fni='out.txt';
- fid=fopen('out.txt','r');
- sf=fscanf(fid,'%f',1);%采样频率
- fi=fscanf(fid,'%f',1);%最小截止频率
- fa=fscanf(fid,'%f',1);%最大截止频率
- t1=fscanf(fid,'%f',1)%时间长度(秒)
- am=fscanf(fid,'%f',1);%最大幅值
- fno=fscanf(fid,'%f',1);%输出数据文件名;
- status=fclose(fid);
- %计算白噪声信号数据长度n
- n=round(sf*t1)+1
- t=0:1/sf:(n-1)/sf;%建立信号离散时间向量t
- nfft=2^nextpow2(n);%大于并接近n的2的幂次方为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,生成幅值谱
- a=2*pi*rand(1,nfft/2);%生成0~2pi的随机数为随机相位
- 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)));%去指定长度的数据并使最大值为指定幅值
- nft=1024;%定义自功率谱的FFt长度
- y1=psd(y,nft,sf,boxcar(nft),nft/2);%计算生成的白噪声信号的自功率谱
- f=0:sf/nft:(nft/2-1)*sf/nft;%定义自功率谱的离散频率向量
- %绘制白噪声信号随时间变化的图形
- subplot(2,1,1);
- plot(t,y);
- xlabel('时间 (s)');
- ylabel(' 幅值 ');
- grid on
- %绘制自功率谱图形
- subplot(2,1,2);
- plot(f(1:nft/4),y1(1:nft/4));
- xlabel('频率 (HZ)');
- ylabel('幅值');
- grid on;
- save force y
- % %打开文件输出生成的白噪声信号数据
- % fid=fopen(fno,'w');
- % for k=1:n
- % fprintf(fid,'%f %f/n',t(k),y(k));%每一行输出两个实型数据,t为时间,y为白噪声信号值
- % end
- % status=fclose(fid);
复制代码 这是书中给出的程序 |