可以进行FFT分析及作图的M文件函数

为了方便进行对已有数据的FFT分析及作图，编制了FFTplot模块，全套模块文件在这里可以下载。它可以对数据型向量(包含时间和幅值)进行FFT分析，并进行绘图。
■ 其中调用了FFTpositive.m和THDcalc.m两个子模块；
■ 读取数据来源可以是1)一般的Scope保存的'Structure with time'类型的变量，也可以是2)Tek示波器直接保存为的CSV文件。选择读取数据来源需要在源代码中作简单的修改。
■ 绘图设置可以在源代码中作简单修改，包括位置、尺寸、标题、轴名称、字体、颜色等。

在这里下载全套模块文件后，可以调用附带的演示数据文件测试效果。

 

验证了FFTplot计算和PowerGUI FFT Analysis的结果是一致的，如下。

以下是其中主模块的代码。

% 20100407              Ver 1.01
% 1.    从Tek示波器保存的csv数据文件中读取数据；
%
% 20100405              Ver 1.00
% 1.    对仿真结果数据进行分析及绘图，绘制信号波形，对信号进行THD分析并作图。
%       数据格式为Scope保存的Structure with time类型
% 2.    调用了FFTpositive.m, THDcalc.m等函数。

clc;
% 初始设置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F1              = 20;                   % 信号基频 !注意和信号的实际基频要设置一致
disFN           = 19;                   % 要显示的最高基波倍次
DATA_FROM       = 1;                    % 0:Workspace; 1:CSV;
CSV_NAME        = 'e:\NOWjobs\ML\Program\DSP\Experiments\ALL0002\F0002CH1.CSV';       % CSV文件名，最好不要有中文
CSV_RANGE       = 'D1..E2500';           % CSV文件数据区域

% 作图设置
AXIS_YMAX         = 0.26%;1;%             % THD波形最大纵坐标值
TITLE_STR        = 'Uan信号THD = ';      % 标题字符串
YLABEL_STR        = '幅值(标幺化)';    %
XLABEL_STR        = '频率(Hz)';
FONTSIZE_AXES   = 12;
FONTSIZE_TITLE     = 14;
FONTSIZE_LABEL     = 12;
BAR_WIDTH          = 0.5;
BAR_COLOR          = 'r';%'b';%

% 提取数据
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if DATA_FROM==0
    % 从Matlab工作空间
    D         = L3;                            % 取工作空间中的数据
    tsim    = D.time;
    X        = D.signals(1,3).values(:,1);
elseif DATA_FROM==1
try
    curpath = cd;
%     fst     = strfind(CSV_NAME,'\');      % 跳转到文首设置的数据文件目录
%     pst     = CSV_NAME(1:(fst(end-1)));
%     pst     = strrep(pst,'\','/');
%     cd(pst);
    [fname,pname]=uigetfile('*.csv','请选择CSV数据文件');    % 通过对话框获得文件路径及名称
    cd(curpath);
    CSV_NAME = [pname,fname];
    D       = dlmread(CSV_NAME,',',CSV_RANGE);      % 直接路径读取csv数据文件
catch
    error('文件读取错误！你是否选择了正确的数据文件路径?');
    return;
end
    tsim    = D(:,1);
    X        = D(:,2);
end

len     = length(tsim);                 % 数据长度
Fs      = 1/(tsim(2)-tsim(1));            % 采样频率

% FFT分析
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 倍频点序号计算
thdFN     = zeros(1,disFN);            % 存放基频和倍频对应的数据点序号
Ft         = Fs/len;                    % 频率间隔
for i=1:disFN
    thdFN(i) = round(F1*i/Ft)+1;    % 基频及倍频点对应的序号
end

% FFT分析
[XP, XF] = FFTpositive(X, Fs);        % 做THD分析,求取频率和幅度
XTHD= THDcalc(XP, XF, F1, disFN);    % 计算THD值
XP    = XP([thdFN]);                    % 取倍频点数据
XF    = XF([thdFN]);
XP    = XP/(XP(1));                    % 幅频以基频幅频值进行标幺

% 作图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 绘图设置位置和尺寸
fig_Width   = 600;
fig_Heigth  = fig_Width*0.6;
scrsz       = get(0,'ScreenSize');  % [left, bottom, width, height]
fig_Position= [(scrsz(3)-fig_Width)/2 ...    % 居中位置
            (scrsz(4)-fig_Heigth)/2 ...
            fig_Width fig_Heigth];
% X轴刻度
xtickstep = 100;                   % X轴坐标刻度单位
xtick0 = 50;
xticknum = ceil(F1*disFN/xtickstep);
xtick = xtick0 : xtickstep : (xticknum*xtickstep);
xticklabel = cell(1,xticknum);
for i=1:xticknum
    xticklabel(1,i) = {num2str((i-1)*xtickstep + xtick0)};
end
% Y轴刻度
ytickstep = 0.02;                   % Y轴坐标刻度单位
yticknum = ceil(AXIS_YMAX/ytickstep);
ytick = 0:ytickstep:AXIS_YMAX;
yticklabel = cell(1,yticknum);
for i=1:yticknum+1
    yticklabel(1,i) = {num2str((i-1)*ytickstep)};
end

% figure('Name', 'Hamonic Spectrum Diagram',...
figure('Name', '谐波频谱图',...
        'MenuBar','figure',...
        'Toolbar','figure',...
        'Color','w',...
        'Position',fig_Position,...
        'NumberTitle','off'...
        );

a1 = axes('FontSize', FONTSIZE_AXES);
h1    = bar(XF, XP, BAR_COLOR);
set(h1,'barwidth',BAR_WIDTH, 'EdgeColor', BAR_COLOR);
% set(get(h1(1),'BaseLine'),'LineWidth',2,'LineStyle',':','Color','red')
title([TITLE_STR, num2str(XTHD*100, '%3.2f'),'% '],'FontWeight','bold','FontSize',FONTSIZE_TITLE);
xlim([0,disFN*F1]);
ylim([0,AXIS_YMAX]);
% set(a1,'XTick',xtick, 'XTickLabel',xticklabel);
set(a1,'YTick',ytick, 'YTickLabel',yticklabel);
set(a1,'XGrid','on')
set(a1,'YGrid','on')
xlabel(XLABEL_STR,'FontSize',FONTSIZE_LABEL);
ylabel(YLABEL_STR,'FontSize',FONTSIZE_LABEL);