基于Matlab与Logistic Regression（逻辑回归）的瓶子密封性检测

DSP大作业内容，做完之后Matlab可能要放一放，把自己一学期收获写下来，留待日后参考。
本设计主要依靠频谱分析和逻辑回归制造简单地二分类器.通过时域阈值将每一次敲击的声音截取，将其转换到频域，取40个500-1500Hz频率采样点作为训练样本，用Matlab中fminunc函数进行训练，得到传输函数Thera，用Thera去判断测试集中的声信号，从而对瓶子密封性做出检测，还有一个比较不错的GUI界面。
首先是对瓶子的声音信号进行截短分段，将每一次敲击的声音提取出来：

global sample_good;
global sample_bad;
global fs;
global theta;
[sounds_good,fs]=audioread('好瓶.wav');%read sound 
[sounds_bad,fs]=audioread('坏瓶.wav');
a1=find(abs(sounds_good)>0.9);%find the louder sound
j=1;
for i=2:length(a1)
  if (a1(i)-a1(j)<1000)
     a1(i)=0;
  else j=i;
  end
end
j=1;
for i=1:length(a1)
  if(a1(i)~=0)
    b1(j)=a1(i);
    j=j+1;
  end
end
a2=find(abs(sounds_bad)>0.9);
j=1;
for i=2:length(a2)
  if (a2(i)-a2(j)<1000)
     a2(i)=0;
  else j=i;
  end
end
j=1;
for i=1:length(a2)
  if(a2(i)~=0)
    b2(j)=a2(i);
    j=j+1;
  end
end
sample_good=zeros((length(b1)-8),2048);%put eight good sample for test
sample_bad=zeros((length(b2)-2),2048);%put two bad samples for test
for i=1:(length(b1)-8)%cut the 2048 points around the louder points as sample
  sample_good(i,:)=(sounds_good((b1(i)-900):(b1(i)+1147)))';
end
for i=1:(length(b2)-2)
  sample_bad(i,:)=(sounds_bad((b2(i)-900):(b2(i)+1147)))';
end

然后将样本通过fft转换到频域并作归一化，取出约500-1500Hz的频率中40个点作为训练样本，样本频率点计算为：原音频采样频率fs=44100Hz，取2048点做fft，取出fft中第30个频率点为：

$$f = \dfrac{44100}{2048}*30 Hz$$

len1=size(sample_good,1);
len2=size(sample_bad,1);
X=zeros(len1+len2,40);
for i=1:len1
    X(i,:)=fft_normalization(sample_good(i,:));
end
for i=1:len2
    X(len1+i,:)=fft_normalization(sample_bad(i,:));
end

function frequency=fft_normalization(time)%时域到频域转换+取645-1485
temp=abs(fft(time));
a1=max(temp);
temp=temp/a1;
frequency=temp(30:69);

将训练样本采到之后，直接调用Matlab的 optimset函数设置参数，调用fminunc函数生成传输函数Thera。以下生成过程数学证明见Courses公开课上机器学习的教学课程（详细）， 逻辑回归模型简介:http://www.cnblogs.com/sparkwen/p/3441197.html?utm_source=tuicool&utm_medium=referral
Matlab逻辑回归简单使用：http://blog.sina.com.cn/s/blog_890c6aa301015mya.html

X=[ones(size(X,1),1) X];
size(X)
y=[ones(len1,1);zeros(len2,1)];
initial_theta = zeros(size(X, 2), 1);
lambda = 0;
options = optimset('GradObj', 'on', 'MaxIter', 300); %number of iteration set as 300
[theta, J, exit_flag] = fminunc(@(t)(costFunctionReg(t, X, y, lambda)), initial_theta, options);

function [J, grad] = costFunctionReg(theta, X, y,lambda) %Regularization coefficient of cost function，
m = length(y);
theta1=theta;
theta1(1)=0;
J=((-y'*log(sigmoid(X*theta))-(1-y)'*log(1-sigmoid(X*theta)))/m)+((lambda*theta1'*theta1)/(2*m));
grad=(X'*(sigmoid(X*theta)-y)/m)+lambda*theta1/m;

function g = sigmoid(z) %S function
g=1./(1+exp(-z));

传输函数Thera创造完成后，用之进行测试。像上述取训练样本同样方法得到测试样本，通过测试样本与传输函数的矩阵相乘，得到二分类的置信度p，通过对p大小的判断就可分辨出该瓶是好瓶还是坏瓶。

global test;
global theta;
test_fre=zeros(1,40);
test_fre=fft_normalization(test);%fft+sampling
fr=35:74;
fr=fr*44100/2048;
axes(handles.iuput);
stem(fr,test_fre,'.');
test_fre=[1;test_fre];
p=sigmoid(test_fre' * theta);%get the confidence interval
if(p>=0.5)
    set(handles.edit1,'String','好瓶')
else
    set(handles.edit1,'String','坏瓶')
end
 set(handles.confidence_level,'String',p)

此程序代码为包含GUI设计代码，无法直接粘贴运行。完整工程文件见链接:wait!
GUI界面：