PSD.rar_matlab psd函数_power spectrum_psd_psd代码 matlab_psd功率谱代码

%% 功率谱绘制 % 所有时间段功率谱 clear clc load('data_sample_2012-07-17 12-VIB.mat'); Data=data(:,1); % Data=xlsread('cj1',1,'a3:z9930'); % Data=xlsread('xiangdui_liu'); % x=detrend(dis1(:,1)); y=detrend(Data(:,1)); % 地震时刻 fs=20; % y=detrend(L(:,1)); nfft=2^nextpow2(length(y)); %index=~isnan(y); %y=y(index); % 其中x是被测信号，window是窗函数，noverlap是相邻两段之间重叠的长度，nfft是FFT变换的长度，Fs是采样频率； % 而pxx是功率谱密度函数，f是频率矢量，w是归一化的角频率矢量。 % 这些参数的默认值：window是汉明窗，信号分为8段，相邻间noverlap是一段长度的一半； % nfft要大于256，或最接近一段长度的2的幂次。 % 当在window参数中没有给出窗函数，只给出一个数值，则用该数值的汉明窗。 % 调用中就是用了L/SmoothingK长的汉明窗；在noverlap中应该是用默认值，即相邻两段间有50%重叠。 % 但如果用L=2^14 ，则分不了段，用pwelch的结果同periodogram的结果相同。 % 在计算自功率谱时是Sx与Sx共轭相乘，得到|Sx|^2，所以为实数，而计算互功率谱时是Sx与Sy共轭相乘，得到的是复数。 window=hamming(nfft/4); noverlap=nfft/8; [pxx, f]=cpsd(y,y,window,noverlap,nfft,fs); % [pxx, f]=cpsd(x,y,window,noverlap,nfft,fs); plot(f,pxx) xlabel('频率 (Hz)','fontsize',20); ylabel('功率谱 ','fontsize',20); set(gca,'fontsize',20) xlim([1 10]); %非地震时刻 %y1=VIB(540000:560000,13); %nfft1=2^nextpow2(length(y1)); %index1=~isnan(y1); %y1=y1(index1); %window1=hamming(nfft1/4); %noverlap1=nfft1/8; %[pxx1, f1]=cpsd(y1,y1,window1,noverlap1,nfft1,fs); %pwelch绘制功率谱 %[pxx1, f1]=pwelch(y1,window1,noverlap1,nfft1,fs); % 非地震时刻 %figure %x_pixels=1400; %y_pixels=900; %figure('units','pixels','position',[0.1*x_pixels,0.1*y_pixels,0.8*x_pixels,0.8*y_pixels]); %set(gcf,'color','w'); %subplot(2,1,1) %plot(f,pxx) %title('北汊桥厦门侧桥塔ACC\_T01\_001','Fontsize',20) %xlabel('频率 (Hz)','fontsize',20); %ylabel('功率谱 (cm^2/s^3)','fontsize',20); %set(gca,'fontsize',20) %xlim([0 10]); %地震时刻前后 %subplot(2,1,2) %plot(f1,pxx1,'color',[0 0 1]) %title('北汊桥厦门侧桥塔ACC\_T01\_001','Fontsize',20) %xlabel('频率 (Hz)','fontsize',20); %ylabel('功率谱 (cm^2/s^3)','fontsize',20); %set(gca,'fontsize',20) %xlim([0 10]); %% 所有时间段功率谱VIC %clc %clear %load('VIC.mat');% % 非地震时刻 %fs=50; %y=VIC(620000:640000,1); %nfft=2^nextpow2(length(y)); %index=~isnan(y); %y=y(index); %window=hamming(nfft/4); %noverlap=nfft/8; %[pxx, f]=cpsd(y,y,window,noverlap,nfft,fs); %地震时刻前后 %y1=VIC(540000:560000,1); %nfft1=2^nextpow2(length(y1)); %index1=~isnan(y1); %y1=y1(index1); %window1=hamming(nfft1/4); %noverlap1=nfft1/8; %[pxx1, f1]=cpsd(y1,y1,window1,noverlap1,nfft1,fs); %%pwelch绘制功率谱 %[pxx1, f1]=pwelch(y1,window1,noverlap1,nfft1,fs); % 非地震时刻 %x_pixels=1400; %y_pixels=900; %figure('units','pixels','position',[0.1*x_pixels,0.1*y_pixels,0.8*x_pixels,0.8*y_pixels]); %set(gcf,'color','w'); %subplot(2,1,1) %plot(f,pxx,'color',[0 0 1]) %title('ACC\_C01\_001\_北汊桥主桥C01垂直拉索方向','fontsize',20) %xlabel('频率 (Hz)','fontsize',20); %ylabel('功率谱 (cm^2/s^3)','fontsize',20); %set(gca,'fontsize',20) %xlim([0 10]); %xlim([0 10]); %地震时刻前后 %subplot(2,1,2) %plot(f1,pxx1,'color',[0 0 1]) %title('ACC\_C01\_001\_北汊桥主桥C01垂直拉索方向','fontsize',20) %xlabel('频率 (Hz)','fontsize',20); %ylabel('功率谱 (cm^2/s^3)','fontsize',20); %set(gca,'fontsize',20) %xlim([0 10]);