【信号处理】一种热电偶信号处理算法

热电偶信号处理，无非两种方法：

数据处理方法
插值法	拟合法
插值表密度决定精度 复杂度集中在搜索插值点	数据密度和拟合方式决定精度 复杂计算耗时大 实现更容易

因为在所设计系统应用中，需要更新插值表，而再导出到PC软件中做拟合，太麻烦（在MCU中实现的拟合也有，例如之前文献学习提到的基于递归B样条最小二乘的方法给出拟合函数，但在单片机里写还是太难了），所以考虑用插值方法。

插值方法时间复杂度主要集中在插值地址搜索，在之前学习的文献中有一种方法：分度表等间隔分布，通过热电势的值简单运算得到查找表地址的思路；这和有序向量插值查找的思路十分相似。

L. Ximin, "A Linear Thermocouple Temperature Meter Based on Inverse Reference Function," 2010 International Conference on Intelligent Computation Technology and Automation, 2010, pp. 138-143, doi: 10.1109/ICICTA.2010.284.

而文献中的这种查找方式有很大的局限性：它必须要求查找表按照Et等间隔、在测量范围内全部有序的给出，这导致，如果想要提高查询表密度以提高精确的，就得全面地提高查找表的密度，占用更大的静态内存，而且需要E-t，t-E两个表。

基于此，我准备利用有序向量插值查找的思路，对这个查找方法进行改写，实现在loglogn的时间复杂度内实现Et到t 的查询、插值，且不必要求查询表等间隔分布，可以在重点测温范围内大密度做表，在其他地方小密度做表:

用MATLAB实现，因为之后可能要配合Simulink做模拟电路、数字数据处理的联合仿真：

先自己写一个插值函数，之后需要改写成C语言：

function [index,Tout] = insert_search(indextable_T,indextable_V,Vin)
    Vsearch=Vin-indextable_V(1);
    id1=Vsearch/(indextable_V(end)-indextable_V(1))*length(indextable_T);
    index=fix(id1);
    if index==0
    stepT=(indextable_T(2)-indextable_T(1));
    stepV=(Vin-indextable_V(1))/(indextable_V(2)-indextable_V(1));
    Tout=indextable_T(1)+stepT*stepV;
    return
    end
    while Vin<indextable_V(index)||Vin>indextable_V(index+1)
        if Vin<indextable_V(index)
        Vsearch=Vin-indextable_V(1);
        id1=Vsearch/(indextable_V(index)-indextable_V(1))*index;
        index=fix(id1);
        elseif Vin>indextable_V(index+1)
        Vsearch=Vin-indextable_V(index);
        id1=Vsearch/(indextable_V(end)-indextable_V(index))*(length(indextable_T)-index);
        index=index+fix(id1);
        end
    end
    stepT=(indextable_T(index+1)-indextable_T(index));
    stepV=(Vin-indextable_V(index))/(indextable_V(index+1)-indextable_V(index));
    Tout=indextable_T(index)+stepT*stepV;

end

这里用到的插值搜索的思想：

 插值得到大概的index，再判断Vin是不是在Vsearch[index，index+1]内，是的话直接线性插值，不是的话再判断是大了还是小了进入[index,end]或者[1,index]递归处理

考虑到所做系统对快速性的需求，所以改写算法的时候把最初的递归改成了while循环迭代，减小了函数栈开销。

测试：

clear; clc; close all

indextable_T = xlsread('K型热电偶分度表.xls',1, 'A3:A143'); %in ℃
indextable_V = xlsread('K型热电偶分度表.xls',1, 'B3:B143'); %分度表数据导入,in mV

V_in=input('模拟热电偶输出电压 -0.392~51.82: ');
T_in=input('模拟冷端补偿温度 -0~1370: ');
[indexT,Tout] = insert_search(indextable_T,indextable_V,V_in)
[indexV,Vout] = insert_search(indextable_V,indextable_T,T_in)

V_tc=(indextable_V(1):0.1:indextable_V(end));
T_tc=(1:length(V_tc));
for i=1:length(V_tc)
    [~,Ttc_now]=insert_search(indextable_T,indextable_V,V_tc(i));
    T_tc(i)=Ttc_now;
end

如果输入定值的V_in或者T_in,可以看到大多数时候2~4次搜索就可以找到地址了，这个效率是非常高的。

再检查输出更高密度V_tc分度表，看上去没什么大问题

 V_tc(256) = 25.5000

 T_tc(256) = 614

原分度表：

600	24.905
610	25.33
620	25.755
630	26.179

---

在考虑另外两个问题：异常值处理和分度表更新

异常值处理会的有以下方法：

3σ检验、Grubbs检验法：  这需要实时计算标准不确定度，之前写在STM32里的不确定度计算函数，由于sqrt函数实现的不够好，总体效果不是很满意，而且单片机里面 u8 char，u32 int  float 和ADC转换的二进制数  相互运算、处理，总不知道哪一步是没转换还是怎么的，例如：

读取的一字节二进制数：

for (i=1;i<=8;i++) 
	{
        j=DS18B20_Read_Bit();
        dat=(j<<7)|(dat>>1);
    }		

    u8 temp;

    TL=Read_Byte(); // LSB   
    TH=Read_Byte(); // MSB  
                        
    tem=TH; //获得高八位
    tem<<=8;    
    tem+=TL;//获得底八位
    tem=(float)tem;//转换

这样从存储器读取一个十六位二进制数，然后转到了 u32 float，32位单精度浮点数。之后如果和double、int之类的做运算，是不是都要先把两边转成float，例如：

float f = 3.4 这样就是double3.4 应该写成： float f = (float)3.4或者写成 float f = 3.4f，总之标准差计算过程还需要进一步完善。

四分位数法：
四分位数一（Q1）数据集前半部分的中位数
四分位数二（Q2）整体集的中位数
四分位数三（Q3）集后半部分的中位数
Q3-Q1计算四分位数Q4范围，Q4乘一个恒定值得到IQR(Interquartile Range)
然后有效数据的下限和上限就是[Q1−IQR,Q3+IQR] 

使用四分位法不用计算标准不确定度，会比较快一些，但是因为有效范围是自定义的，就是说不想σ或者Grubbs法这样确定异常值的把握有多少（99.75%，比如）

---

分度表更新，有两个问题：

一是选用何种数据结构存储分度表，二是如何改变本地存储，就是变了之后单片机reset一下数据不会复原。存flash感觉不太安全，所以存了sd卡。

至于分度表的数据结构，简单的两个数组我觉得是不合适的（虽然这最简单而且常用），考虑到分度表更新，数组的话更新限制太大，在一个位置插入T，就要在另一数组相同索引加入E，插入过程也就是后面挪一下，但是数据安全完全没有保障，就是说如果在更新的时候出了一些意外，比如insert(T,E）,T不小心输大十了十倍，存最后去了，那E是跟着T的索引还是自己再寻找插值地址?

如果是E跟着T索引，那要改正只能说在数组里进行delete(T)，删除，后面全部前移，再delete（E）,如果不想这么复杂，在插入的时候得加一个判断，判断E跟着T的插值地址，其插值表两边的值是否合适。 这使得更新的逻辑特别复杂。如果是E，T分布自己寻找插值地址更新，那就失去了单调检验的功能：输入的（E，10T）分别插值到不同的索引，那么也是得定向delete刚插入的，没有发现插入错误的话就得重新制表了，这非常的不合理。

更新困难只是一小部分，另一重要原因是：如果存储的数组被不小心错误更改（因为要提供更新表的功能所以肯定要做修改的接口），例如 E T 两个数组删了索引没对上的一组元素，那么整么表在某个索引后就是错误的了，而且这种情况也能通过顺序排列检查！

概括：把分布表做为两个数组结构不适合更新，其不仅跟新复杂，而且安全性没有保障。要保障安全、可靠，需要非常非常复杂的逻辑判断防止错误操作，注意事项太多，很容易出问题。

改进的方案，一是做一个struct结构、Binary Search Tree结构或者 Class:ThermocoupleData ，每个元素(节点)存两个float数据和指针。 如果是struct结构可以用指针数组指向n个struct，完成分度表建立；如果是二叉搜索树，直接用搜索树Class模板，然后改一下树节点状态即可。

用指针数组的好处是分度表在struct内严格一一对应，这样数据安全有了保障，无论如何不会出现错误对应的数据。而插入、删除、搜索的复杂度都是和数组类似的；用二叉树结构存储，看起来更厉害一些，而且插入、删除、搜索都可以在log(n)复杂度内实现，但是上面写了半天的向量插值搜索就白写了。而且一个麻烦是Class 是C++适应，转成C头文件也挺累的。

所以目前考虑的方案：用struct+指针vector 存储分度表。