psd.m

自己编写的PSD批处理程序

function psd(fs,colPos,startpos,compute_length)
%
%功率谱密度算法,    Write by Zheng guibo   08-05-11

%fs               : sample rate
%startpos         : start position of time series
%compute_length   : number of time series
%colPos           : data column 

%-----------------inital program--------------------
show_f       = 40;   %控制显示频率段
file_ext     = '.xls';           %要计算的文件类型,根据数据文件修改为'.xls'或'.txt'
%%*****注意，列数是指有效数据列，时间列不计算在内，
%%******如果有时间列，则第二列为本程序所指第一列
const_colPos  = 4;            %%数据列位置
const_fs     = 400;          %%默认采样频率
const_start  = 4000;         %%默认采样起点
const_length = 4096;         %%默认计算长度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nargin == 0),
    fs        = const_fs;
    colPos    = const_colPos;
    startpos  = const_start;
    compute_length= const_length;    
    %error('At least sample rate (fs) required');
end;

if (nargin == 1),
    colPos   = const_colPos;
    startpos = const_start;
    compute_length= const_length;    
end;
if (nargin == 2),
    startpos = const_start;
    compute_length= const_length;    
end;
if (nargin == 3),
    compute_length= const_length;    
end;


%%指定一个不存在的文件可以使默认路径为桌面
sel_directory_name = uigetdir('E:\zheng\2008.4.12倾斜数据\751度','请选择数据文件所在路径');
d = dir(strcat(sel_directory_name,strcat('\*', file_ext)));
ac_files = 0;
%d = dir('*.xls');
str = {d.name};
if (isempty(str)==1),
    warndlg('选择目录没有数据文件，请重新选择！','目录选择有误'); 
    return;
else  
    for i = 1:size(str,2)

        filename = str(i);

        k = findstr(char(filename), 'AC');     %只计算AC数据，如果同时计算DC，注释本行，并注释if语句

        if (isempty(k)==0),
            ac_files = ac_files+1;
            fullname = strcat(sel_directory_name,'\',str(i));  

            A = xlsread(char(fullname));              
            %如果不指定计算长度，则默认计算长度为从起点到最后
            if (nargin < 4),    
                compute_length = length(A);
                colData = A(startpos + 1: compute_length, colPos);                %取出到最后一行所有数据
            else
                colData = A(startpos + 1: startpos + compute_length , colPos);    %取出指定长度数据
            end;
%             nfft=1024; 
%             window=boxcar(100); %矩形窗 
%             window1=hamming(100); %海明窗 
%             window2=blackman(100); %blackman窗 
%             noverlap=20; %数据无重叠 
%             range='half' %频率间隔为[0 Fs/2]，只计算一半的频率 
            [pp,ff] = pwelch(colData,[],[],[],fs);

            h=figure(1);
            
            %显示40Hz以下部分

            index = find(ff>=show_f); 
            plot(ff(1:index(1)),(pp(1:index(1))));%10*log10
            set(gca,'Xlim',[0,show_f]);

            xlabel({'Frequence(Hz)';char(filename)},'Interpreter','none');   
            ylabel('Power(db)');
            title('PSD');          
            %找出最大峰值对应频率
            [aa,p_index]=max(pp);
            ss=sprintf('Peak_Frequence= %0.2f Hz',ff(p_index(1)));            
            text('Position',[show_f/2,aa*0.9],'Interpreter','none','string',ss,'VerticalAlignment','bottom');
            
            if ~exist(strcat(sel_directory_name,'\PSD result')) 
                mkdir(strcat(sel_directory_name,'\PSD result')); 
            end
            paths = [sel_directory_name,'\PSD result\'];
            savename  = strrep(filename, file_ext, '.jpg');
            saveas(gcf,strcat(paths,char(savename)));
            close;
            clear pp;
            clear ff;
        end;
    end;
end;

warndlg('计算完成');