在使用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 (option 大专栏 Matlab - 常用函数集锦al)
%
% 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;