Matlab - 常用函数集锦

在使用matlab进行信号处理和图形绘制过程中，某些函数被频繁调用，所以有必要将这些常用函数进行总结归类。

滤波函数

低通滤波

function [filtered_signal,filtb,filta]=lopass_butterworth(inputsignal,cutoff_freq,Fs,order)
% Low-pass Butterworth filter
% [filtered_signal,filtb,filta] = lopass_butterworth(inputsignal,cutoff_freq,Fs,order)
%
% This is simply a set of built-in Matlab functions, repackaged for ease of
% use by Chad Greene, October 2012.
%
% INPUTS:
% inputsignal = input time series
% cutoff_freq = filter corner frequency
% Fs = data sampling frequency
% order = order of Butterworth filter
%  
% OUTPUTS:
% filtered_signal = the filtered time series
% filtb, filta = filter numerator and denominator (optional)
%
% EXAMPLE 1:
% load train
% t = (1:length(y))/Fs;
% y_filt = lopass_butterworth(y,900,Fs,4); % cut off at 900 Hz
% figure
% plot(t,y,'b',t,y_filt,'r')
% xlabel('time in seconds')
% box off
% legend('unfiltered','filtered')
% sound(y,Fs)      % play original time series
% pause(2)         % pause two seconds
% sound(y_filt,Fs) % play filtered time series

nyquist_freq = Fs/2;  % Nyquist frequency
Wn=cutoff_freq/nyquist_freq;    % non-dimensional frequency
[filtb,filta]=butter(order,Wn,'low'); % construct the filter
filtered_signal=filtfilt(filtb,filta,inputsignal); % filter the data with zero phase

高通滤波

function [filtered_signal,filtb,filta]=hipass_butterworth(inputsignal,cutoff_freq,Fs,order)
% High-pass Butterworth filter
% [filtered_signal,filtb,filta] = hipass_butterworth(inputsignal,cutoff_freq,Fs,order)
%
% This is simply a set of built-in Matlab functions, repackaged for ease of
% use by Chad Greene, October 2012.
%
% INPUTS:
% inputsignal = input time series
% cutoff_freq = filter corner frequency
% Fs = data sampling frequency
% order = order of Butterworth filter
%
% OUTPUTS:
% filtered_signal = the filtered time series
% filtb, filta = filter numerator and denominator (optional)
%
% EXAMPLE 1:
% load train
% t = (1:length(y))/Fs;
% y_filt = hipass_butterworth(y,900,Fs,4); % cut off at 900 Hz
% figure
% plot(t,y,'b',t,y_filt,'r')
% xlabel('time in seconds')
% box off
% legend('unfiltered','filtered')
% sound(y,Fs)      % play original time series
% pause(2)         % pause two seconds
% sound(y_filt,Fs) % play filtered time series

nyquist_freq = Fs/2;  % Nyquist frequency
Wn=cutoff_freq/nyquist_freq;    % non-dimensional frequency
[filtb,filta]=butter(order,Wn,'high'); % construct the filter
filtered_signal=filtfilt(filtb,filta,inputsignal); % filter the data with zero phase

带通滤波

function [filtered_signal,filtb,filta]=bandpass_butterworth(inputsignal,cutoff_freqs,Fs,order)
% Bandpass Butterworth filter
% [filtered_signal,filtb,filta] = bandpass_butterworth(inputsignal,cutoff_freq,Fs,order)
%
% This is simply a set of built-in Matlab functions, repackaged for ease of
% use by Chad Greene, October 2012.
%
% INPUTS:
% inputsignal = input time series
% cutoff_freqs = filter corner frequencies in the form [f1 f2]
% Fs = data sampling frequency
% order = order of Butterworth filter
%
% OUTPUTS:
% filtered_signal = the filtered time series
% filtb, filta = filter numerator and denominator (optional)
%
% EXAMPLE 1:
% load train
% t = (1:length(y))/Fs;
% y_filt = bandpass_butterworth(y,[800 1000],Fs,4); % cut off below 800 Hz and above 1000 Hz
%
% figure
% plot(t,y,'b',t,y_filt,'r')
% xlabel('time in seconds')
% box off
% legend('unfiltered','filtered')
% sound(y,Fs)      % play original time series
% pause(2)         % pause two seconds
% sound(y_filt,Fs) % play filtered time series

nyquist_freq = Fs/2;  % Nyquist frequency
Wn=cutoff_freqs/nyquist_freq;    % non-dimensional frequency
[filtb,filta]=butter(order,Wn,'bandpass'); % construct the filter
filtered_signal=filtfilt(filtb,filta,inputsignal); % filter the data with zero phase

带阻滤波

function [filtered_signal,filtb,filta]=bandstop_butterworth(inputsignal,cutoff_freqs,Fs,order)
% Band-stop Butterworth filter
% [filtered_signal,filtb,filta] = bandstop_butterworth(inputsignal,cutoff_freqs,Fs,order)
%
% This is simply a set of built-in Matlab functions, repackaged for ease of
% use by Chad Greene, October 2012.
%
% INPUTS:
% inputsignal = input time series
% cutoff_freqs = filter corner frequencies in the form [f1 f2]
% Fs = data sampling frequency
% order = order of Butterworth filter
%
% OUTPUTS:
% filtered_signal = the filtered time series
% filtb, filta = filter numerator and denominator (optional)
%
% EXAMPLE 1:
% load train
% t = (1:length(y))/Fs;
% y_filt = bandstop_butterworth(y,[800 1000],Fs,4); % cut off below 800 Hz and above 1000 Hz
%
% figure
% plot(t,y,'b',t,y_filt,'r')
% xlabel('time in seconds')
% box off
% legend('unfiltered','filtered')
% sound(y,Fs)      % play original time series
% pause(2)         % pause two seconds
% sound(y_filt,Fs) % play filtered time series

nyquist_freq = Fs/2;  % Nyquist frequency
Wn=cutoff_freqs/nyquist_freq;    % non-dimensional frequency
[filtb,filta]=butter(order,Wn,'stop'); % construct the filter
filtered_signal=filtfilt(filtb,filta,inputsignal); % filter the data with zero phase

绘图函数

function [  ] = setPlot( varargin )
% setPlot()
% setPlot(title)
% setPlot(title,xlabel)
% setPlot(title,xlable,ylabel)
% setPlot(title,xlabel,ylabel,xlim)
% setPlot(title,xlable,ylabel,xlim,ylim)

narginchk(0,5);     % 判断输入参数是否足够

grid on;
axis tight;

if nargin>=1
    title(varargin{1});
end

if nargin>=2
    xlabel(varargin{2});
end

if nargin>=3
    ylabel(varargin{3});
end

if nargin>=4
    xlim(varargin{4});
end

if nargin>=5
    ylim(varargin{5});
end

end

信号处理函数

频谱分析

function [freq,amp]=fft_signal(signal,fs,N)
% Spectrum analysis
% INPUTS:
% signal = input time series
% fs = data sampling frequency
% N = data length of signal
%
% OUTPUTS:
% freq = frequency of Spectrum
% amp = amplitude of Spectrum
%
% EXAMPLE 1:
% fs = 100;
% N = fs*10;
% t = (0:N-1)/fs;
% y = sin(2*pi*10*t);
% [freq,amp] = fft_signal(y,fs,N);
% plot(freq,amp);

amp = 2*abs(fft(signal))/N;        % 求取信号的幅度谱
amp = amp(1:fix(length(amp)/2));       % 截取有效部分
freq=(0:length(amp)-1)*fs/N;       % 横坐标代表频率
end

幅值分布

function [ amp,dist ] = ampDist( signal,sectionNum )
% Calculate the amplitude distribution of the signal
% INPUTS:
% signal : The signal to be analyzed
% sectionNum : Number of segments
%
% OUTPUTS:
% amp : Amplitude after segmentation
% dist :Amplitude distribution
%
% EXAMPLE 1:
% fs = 1000;
% N = fs*100;
% y = wgn(1,N,10);     % 高斯白噪声
% [amp,dist] = ampDist(y,500);
% bar(amp,dist);

yMin = min(signal);
yMax = max(signal);

amp = linspace(yMin,yMax,sectionNum);
dist = hist(signal,amp);
dist = dist./length(signal);
end

LMS最小均方算法

function [ y_error, y_filter ] = LMS( x_input,x_dest,M,u )
% LMS 最小均方算法
% INPUTS：
% x_input    原始信号
% x_dest     期望信号
% M          阶次
% u          步长因子
%
% OUTPUTS：
% y_error    误差信号
% y_filter   滤波器信号输出
%
% EXAMPLE 1:
% load train
% t = (1:length(y))/Fs;
% M = 2; u = 0.5;
% y_dest = (max(y)-min(y))/2*cos(2*pi*18*t);      % 参考信号
% [y_error,y_filter] = LMS(y,y_dest,M,u);
% plot(t,y_error,t,y_filter);

N = length(x_input);
y_filter = zeros(1,N);
y_error  = zeros(1,N);
h = zeros(1,M);

for k=M:N
    h_old = h;
    y_filter(k) = x_dest(k:-1:k-M+1)*h_old';
    y_error(k)  = x_input(k) - y_filter(k);
    h = h_old + 2*u*y_error(k)*x_dest(k:-1:k-M+1);
end

end

EMD经验模态分解

function imf = emd(x)
% Empiricial Mode Decomposition (Hilbert-Huang Transform)
% imf = emd(x)
% Funcs : ismonotonic, isimf, getspline, findpeaks

x   = transpose(x(:));      % 将x变为一维向量
imf = [];
while ~ismonotonic(x)
   x1 = x;
   sd = Inf;
   cnt=0;
   while (sd > 0.1) || ~isimf(x1)
      s1 = getspline(x1);
      s2 = -getspline(-x1);
      x2 = x1-(s1+s2)/2;

      sd = sum((x1-x2).^2)/sum(x1.^2);
      x1 = x2;
      cnt=cnt+1;
   end
%    cnt

   imf{end+1} = x1;
   x          = x-x1;
end
imf{end+1} = x;

% FUNCTIONS
% 判断信号的单调性
function u = ismonotonic(x)

u1 = length(findpeaks(x))*length(findpeaks(-x));
if u1 > 0
    u = 0;
else
    u = 1;
end

% 判断信号是否满足IMF条件
% 条件：极大值点数和极小值点数之和与过零点数相等或相差1？
function u = isimf(x)

N  = length(x);
u1 = sum(x(1:N-1).*x(2:N) < 0);
u2 = length(findpeaks(x))+length(findpeaks(-x));
if abs(u1-u2) > 1
    u = 0;
else
    u = 1;
end

% 使用三次样条函数，得到包络线
function s = getspline(x)

N = length(x);
p = findpeaks(x);
s = spline([0 p N+1],[0 x(p) 0],1:N);

% 寻找极大值点
function n = findpeaks(x)
% Find peaks.
% n = findpeaks(x)

n    = find(diff(diff(x) > 0) < 0);
u    = find(x(n+1) > x(n));
n(u) = n(u)+1;