RL反卷积都干了些啥

RL反卷积是通过迭代的方法来优化损失函数,利用噪声性质和成像系统的点分布函数(Point Spread Function)对图像进行去噪的一种方法。 损失函数是啥样的呢?

且有:

  • ,损失函数,为最小化的目标
  • ,不受椒盐噪声和点分布函数模糊效应影响的理想原始图像
  • ,Point Spread Function,点分布函数
  • ,实际拍摄到的图像
  • ,二维卷积
  • ,按元素相乘

既然要优化,那么就对O求导:

利用迭代法,就可以使代价函数最小,如下:

简化可得:

按照这个迭代式,一般RL反卷积算法迭代个十几次就差不多了。

但是!

为什么要这么干呢?代价函数是怎么来的呢? 这和光学设备的成像的噪声性质有关系。

数码噪点的泊松分布

由于光线具有波粒二象性,在进行数字成像时,在像素点上的光子数量在概率上呈泊松分布。

对于某个像素点而言,曝光时的光电子数量满足泊松分布:

其中,表示像素接收到i个光电子的概率,而表示像素接收到的光电子数量的期望值(平均值)。

问题的转化:极大似然估计

考虑整个图像的联合概率为每个像素的概率的乘积:

且有:

  • :整个图像每个像素的联合概率
  • :实际拍摄到的图像数据,包含噪声,M行N列
  • :无泊松噪声的理想图像数据,M行N列

那么很自然的,才用极大似然估计,对上式取对数,化乘为加:

其中,是常数项,而对于,有:

去掉到最大似然估计函数前面的负号以及常数项,代入上式就可以得到最小化目标的损失函数:

小结

最大似然估计大法好。一些机器学习方法里面的贝叶斯方法也用到了最大似然估计来求参,真是有种触类旁通的感觉呢。

参考文献

  1. Elements of Richardson-Lucy Deconvolution
  2. The Richardson-Lucy Algorithm
  3. Photographic Sensor Simulation
  4. 入门——理解数码噪点:产生原因及对拍摄的指导意义
  5. Richardson W H. Bayesian-based iterative method of image restoration[J]. JoSA, 1972, 62(1): 55-59.

unmet dependency

本着能省一事就省一事的原则,我装了Ubuntu 20.04 LTS,可是在安装ROS noetic的时候还是出现了unmet dependency的问题。 默认的--fix-missing是没有用的,看看哪些包出问题了就用下面的命令强制覆写一下就OK了。 这些包一般会缓存在/var目录下的子目录里面,报错信息里面会有说明。

1
2
sudo dpkg -i --force-overwrite <filename>
sudo apt install --fix-missing

ros_readbagfile.py

这个文件呢本来是为了更有效率地阅读bag文件的,据wiki说比rostopic要快7倍。 在实际使用过程中出现了依赖问题,提示缺少动态链接库之类的。 从python2转换到python3解决了这个问题。

1
2
# 改下文件头的解释器设置
#!/usr/bin/python3
1
2to3 -f has_key -w ros_readbagfile.py

network configuration

在测试多个笔记本的通信时,发现master已经publish消息了,但是另外一台笔记本却没有收到。这是因为localhost对应的IP地址为127.0.0.1,只能在本机访问。设置master笔记本的ROS_IP环境变量即可。

1
export ROS_IP='192.168.***.***'

roswtf: what the fuck?

在出现问题的时候,用这个命令排查一下,还是挺方便的:)

Reference

  1. ROS: Troubleshooting
  2. ROS: NetworkSetup

生活意义

我本来是不觉得生活本身有什么意义。看完Gorden Mathews的前两集讲座之后,也依然觉得生活本身并没有什么意义。不过确实认同工作最终是不适合成为生活的意义的。

如果大多数人追寻生活的意义,那么社会就会瓦解。现实是大多数往往都会自动适应社会的惯例,自主地或随波逐流地将自己代入社会分工中去,并实行自我欺骗赋予生活以意义。

从工作中获取地意义往往在于两点:

  1. 可衡量的进步:职称、薪资、销售额等等
  2. 创造

创造性的工作总是稀缺的,工作内容中的创造性部分也是稀缺的。而可衡量的进步往往更加实在、具体。

很多名言、书籍都会告诉人们应该如何生活,追寻怎样的人生,泛滥的引用充斥在各种角落。我只想引用一句散落在偏僻角落的话来作结:

当我还是一个相当早熟的少年的时候,我就已经深切地意识到,大多数人毫无休止地追逐的那些希望和努力都是毫无价值的。而且我不久就发现了这种追逐的残酷,这在当年较之今天是更加精心地用伪善和漂亮的字句掩饰着的。

1. 简介

本文档主要根据LeMaRiva的博客内容整理并基于最新的rt分支代码复现,简述以下2个方面:

  1. 如何为树莓派编译实时内核
  2. 测试实时内核的进程调度延时性能

实验所需硬件:

  1. 树莓派3B
  2. 8G TF卡

2. 编译

1.获取代码和工具:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
# From LeMaRiva
mkdir rpi-rt
cd rpi-rt
mkdir rt-kernel # 用于存放编译产生的文件
cd rt-kernel
mkdir boot
cd ..
# 获取树梅派源码
git clone https://github.com/raspberrypi/linux.git
# 获取编译工具
git clone https://github.com/raspberrypi/tools.git
cd linux
git checkout rpi-4.14.y-rt # 截至2019年6月16日最新的rt分支

# 配置变量
export ARCH=arm
export CROSS_COMPILE=~/rpi-rt/tools/arm-bcm2708/gcc-linaro-arm-linux-gnueabihf-raspbian/bin/arm-linux-gnueabihf-
export INSTALL_MOD_PATH=~/rpi-rt/rt-kernel
export INSTALL_DTBS_PATH=~/rpi-rt/rt-kernel

# 对RPi 2/3 B以及RPi Compute Module的配置(Compute Module本人测试过,使用该配置)
export KERNEL=kernel7
make bcm2709-defconfig
# 对RPi A(+), B(+), Zero
export KERNEL=kernel
make bcmrpi_defconfig

vim .config
# 确认一下字段的配置
CONFIG_HIGH_RES_TIMERS=y # 默认已配置为y
CONFIG_PREEMPT_RT_FULL=y # 默认已配置为y
CONFIG_HZ_1000=y # 需要手动更改
CONFIG_HZ=1000 # 需要手动更改

2.编译内核

1
2
3
4
5
make -j8 zImage
make -j8 modules
make -j8 dtbs
make -j8 modules_install
make -j8 dtbs_install

编译完成后,rt-kernel文件夹中结构为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
.
├── boot
├── lib
├── overlays
├── bcm2708-rpi-0-w.dtb
├── bcm2708-rpi-b.dtb
├── bcm2708-rpi-b-plus.dtb
├── bcm2708-rpi-cm.dtb
├── bcm2709-rpi-2-b.dtb
├── bcm2710-rpi-3-b.dtb
├── bcm2710-rpi-3-b-plus.dtb
├── bcm2710-rpi-cm3.dtb
├── bcm2835-rpi-a.dtb
├── bcm2835-rpi-a-plus.dtb
├── bcm2835-rpi-b.dtb
├── bcm2835-rpi-b-plus.dtb
├── bcm2835-rpi-b-rev2.dtb
├── bcm2835-rpi-zero.dtb
├── bcm2835-rpi-zero-w.dtb
├── bcm2836-rpi-2-b.dtb
└── bcm2837-rpi-3-b.dtb

注意:此处可能产生警告,有关编译时间对应关系的,可以忽略。

1
2
3
4
./scripts/mkknlimg ./arch/arm/boot/zImage $INSTALL_MOD_PATH/boot/$KERNEL.img
cd $INSTALL_MOD_PATH
tar czf ../rt-kernel.tgz *
scp rt-kernel.tgz pi@<IPAddress>:/tmp

3.替换内核 进入树梅派终端:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
cd /tmp
tar xzf rt-kernel.tgz
cd boot
sudo cp -rd * /boot/
cd ../lib
sudo cp -rd * /lib/
~$ sudo nano /boot/config.txt
## 树梅派 3B 或者 Compute Module
kernel=vmlinuz-4.14.91-rt49-v7+
# 树梅派A(+), B(+), Zero
kernel=vmlinuz-4.14.91-rt49+

sudo reboot
```

如何查看编译内核的版本呢?进入rpi-rt/rt-kernel/lib/modules文件夹下,即可通过文件夹名称看到所编译内核的具体版本号。 [LeMaRiva博客](https://lemariva.com/blog/2018/02/raspberry-pi-rt-preempt-tutorial-for-kernel-4-14-y)中提到了需要在cmdline.txt中进行额外的设置。不过我在3B和Compute Module中都未添加那些额外的配置,测试未发现问题,可能是因为在最新的实时内核中已经修复了bug。重启后查看内核版本号即可确认内核替换成功:

```bash
uname -r
4.14.91-rt49-v7+

3. 验证与测试

可以采用rt-tests工具中的cyclictest对于内核的进程响应延时进行测试。参考Tiejun Chen的分享的脚本,可以自动测试并调用gnuplot生成测试结果图像。

下载编译rt-tests:

1
2
3
4
5
6
7
8
9
10
11
12
13
sudo apt-get install build-essential # libnuma-dev在raspbian stretch中找不到,也不需要
git clone git://git.kernel.org/pub/scm/utils/rt-tests/rt-tests.git
cd rt-tests
git checkout stable/v1.0
make all
sudo make install
```

自动化测试脚本:

```bash
wget https://www.osadl.org/uploads/media/mklatencyplot.bash
sudo apt install gnuplot # 该脚本使用gnuplot

为了节约测试验证的时间,我把该脚本中的次数由1e8改成了1e7。这可能导致样本不够,但是也足够表现实时内核的性能。 实测在默认内核下,最大延时为422us;而在实时内核下,最大延时为118us。

Latency of Raspbian Latency of RT kernel

4. 参考

  1. LeMaRiva's Blog
  2. Tiejun Chen 在OpenIoT Summit of North America上的分享

在当初学《自动控制原理》的时候,并没有去理解公式背后的意义,囫囵吞枣式地做题、计算、画图。现在想来,传递函数的计算、根轨迹的绘制甚至是伯德图的快速绘制方法都记不得了,但是在实践和仿真中,对于概念的理解更深入了些,所以写下。本文从伯德图说起,讲了以下3个方面:

  1. 伯德图表示什么
  2. 开环伯德图表示什么
  3. 伯德积分定理意味着什么

伯德图表示什么

伯德图是从频域的概念来描述系统的性质、表现的。如果我们考虑最简单的系统——比例环节:

时,其传递函数如下图所示:

比例环节伯德图

那么这张图代表什么呢?代表着任意频率的正弦信号通过比例环节之后输出的信号的幅值依然不变——增益为0dB;同时输出信号和输入信号是同相位的,没有相位延迟。

那么如果我们考虑常见的一阶惯性环节呢:

时,一节惯性环节的伯德图是什么样子的呢?

一阶惯性环节伯德图

可见,一阶惯性环节在低频处的增益近似为0dB,而随着频率的增加,增益逐渐下降,而相位也随着频率的增加而滞后,相位滞后的极限为90°。其实一阶惯性环节这样的传递特性对应于电阻电容串联的电路,系统输入为施加在串联阻容上的电压;系统输出为电容两端的电压。如果输入电压为直流电压,那么电容的稳态电压就是直流电压;而如果施加低频正弦电压,那么输出的信号的衰减和相位滞后就会很小;而如果施加高频正弦信号的话,那么输出信号的幅度相对于输入信号就会产生衰减,相位滞后也会随之增大。

那么衰减多少呢?滞后多少呢?对于频率为的正弦信号,其增益和相位滞后如下式所示。

开环伯德图——控制系统稳定分析

考虑下图显示的典型的反馈控制系统框图,如果将反馈处切开,那么从给定(Setpoint 或者 Reference)到反馈的传递函数实际上就是自控概念中的开环传递函数了。

闭环控制系统

因为线性定常系统的性质——叠加性+齐次性——不同频率的输入信号的响应也可以叠加。所以很多时候为什么我们看伯德图呢?因为伯德图反映了系统输入什么样的正弦波之后,会出来什么正弦波呀。而输入信号一般都是可以分解为不同频率和相位的正弦波的叠加啊(想想傅里叶变换?)。

结合上图,对于正弦信号而言,开环传递函数的幅频特性实际上反映的是反馈回来的信号的放大倍数(也可能是缩小)和相角延迟,这个和系统的稳定性有着密切的关系。如果反馈回来的信号为-180°,并且增益大于等于1,那么系统就自然不会稳定。然而在实际系统中,由于系统功率(能量)或者非线性环节的限制、噪声的影响等等,系统的状态可能表现为自激振荡。

而这就是对数稳定判据的一种理解方式了,其实对数稳定判据和奈奎斯特稳定判据殊途同归。

水床效应(waterbed effect)

水床效应实际上就是考虑输出扰动的情况下,系统对于特定频率扰动抑制的局限——如果某种频率范围的扰动被抑制了,那么必然有其他频率范围的扰动被增大了。

为什么呢?伯德敏感度积分如下所示:

其中, 是系统的环路增益,p_{k}是环路增益L(s)的不稳定极点。

如果

  1. 没有不稳定的极点(确保上式被减数为0)
  2. 极点比零点多至少2个(确保上式被减数为0)

那么,就有:

所以说,对于扰动的闭环传递函数而言,某些频率范围的增益的抑制(减小),必然会导致其他某些频率范围增益的增大。这里按下去,那里就会鼓起来,即水床效应咯。

毕竟根据“天下没有免费的午餐”定理,获得任何效用都是要付出额外的代价的:)

参考

  1. Wikipedia: Bode plot
  2. Wikipedia: Bode's sensitivity integral

1.测试环境

  • 操作系统:Windows 10, 64bit
  • 编译器:MinGW64
    • 32位系统请下载对应的32位版本
  • MATLAB2017a
    • 2017a之后会略有不同

2. 配置环境与流程

  1. 下载编译器MinGW64,对于64位系统而言,需要下载个安装器,可以安装相应版本的gcc编译器,我选择的是gcc-4.9.4,在测试环境中可以正常运行

    1. 安装器是mingw-w64-install,安装时注意选择64位,同时安装完了还有一个坑:空格。
    2. 下载器默认的安装路径"C:Files-w64_64-4.9.4-win32-seh-rt_v5-rev0"是包含有空格的,MATLAB并不能识别该路径,所以可以把mingw64移动到不含空格的路径下,比如"C:"。
  2. 验证一下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
setenv('MW_MINGW64_LOC','C:\mingw64') % 先配置环境
mex -setup
```

输出如下图所示,这样MinGW编译器就配置好了。

![配置成功输出](config_success.jpg)


## 3. a + b 编译测试

```MATLAB
setenv('MW_MINGW64_LOC','C:\mingw64')

myFunctionName = 'SFunctionName'

def = legacy_code('initialize');
def.SFunctionName = myFunctionName ;
% 注意,输入和输出的变量的名称只能是u1,u2...和y1,y2...
def.OutputFcnSpec = 'double y1 = add_2(double u1, double u2)';
% 设置用到的头文件或者源文件
def.HeaderFiles = {'add_2.h'};
def.SourceFiles = {'add_2.c}
def.InitializeConditionsFcnSpec = 'buffer_init()';

legacy_code('sfcn_cmex_generate', def);
legacy_code('compile', def) ;

以上代码可以在MATLAB中快捷编译,并且生成Simulink可用的文件"add_2.mex64",在Simulink中加入S-Function模块,并且将名字设置为"add_2"即可。

add_2.c就可以直接写一个两个数相加的程序验证一下就好了; add_2.h写个函数声明也就可以了。

1
2
3
4
5
6
7
8
9
10
11
12
// add_2.c
#include "add_2.h"

double add_2(double a, double b){
return a + b ;
}

// add_2.h
#ifndef ADD_2_H
#define ADD_2_H
double add_2(double a, double b);
#endif

4. 参考

  1. Mathworks关于空格的说明
  2. MATLAB 2017后版本注意事项
  3. 印度小哥的解说

1. Introduction to Proportional-Resonant Controller

First, let's see the transfer function of PR controller and its bode diagram:

where the is the resonant frequency.

Bode diagram of PR Controller

As can be seen in the bode diagram, the gain at 50Hz (314 rad/s) is infinite, so PR controller can be used to track reference of specific frequency. However, the reference frequency is not always a constant, e.g. the frequency of electricity grid. In practice, Quasi-Proportional-Resonant(QPR) can be used to solve this problem.

The transfer function of QPR and its bode diagram:

Bode diagram of QPR Controller

where the addition of reduces the gain at resonant frequency but increase the band width around resonant frequency.

You can change the parameter and run the demo code below to see how it changes the shape of bode diagram.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
s = tf([1, 0], 1) ;

Kp = 15 ;
Kr = 2000 ;
wr = 50 * 2 * pi ;
wi = pi ;

G_PR = Kp + Kr * s / (s^2 + wr^2)
figure
bode(G_PR)

G_QPR = Kp + 2 * wi * Kr * s / (s^2 + 2 * wi * s + wr^2)
figure
bode(G_QPR)

2. Implementation in C

The code is in the repository Controller, feel free to use it for your own application.

The structure of QPR controller implemented is as below:

PR Implementation

As shown in this paper, the structure of QPR controller can be implemented using 2 integrators.

The implemented structure above has 2 advantages:

  • Easy to implemented, only 2 integrators used.
  • Autonomous resonant frequency adaption: can be fed into the controller and real-time resonant frequency adaption is achieved.

You can run this demo to see the effectiveness of this structure

  • is set to be 0 and to be 10
  • is set to be

The output response is as below:

Simulation Result

As you can see in the result above, the input signal is a sine wave with an amplitude of 1 and the output with an amplitude of 10, corresponding to the parameter-- is 0 and is 10.

3. Reference

  1. Effects of Discretization Methods on the Performance of Resonant Controllers
  2. 比例谐振控制算法分析

离散化方法的笔记

最近看论文中发现,分析系统的稳定性问题总免不了在s域先分析一下,然后再在z域上分析一下。可是很多论文对于离散化方法却略去不提,直接给出推导出的结果。所以笔者就寻找总结了一下各种离散化方法,总结一下。学通信或者测量的同学,专门学过数字信号处理的话,应该很熟悉的。

1. 后向差分

后向差分的映射关系为:

从s平面到z平面上来看的话,这实际上是将s平面的左半平面映射到z平面的以(0.5,0)为圆心,半径为0.5的圆中。所以原来稳定的系统,经过后向差分之后,一定是稳定的,因为该圆在单位圆内;而原来不稳定的系统,有可能被映射到单位圆中,即可能存在原来不稳定的系统经过后向差分离散化之后,是稳定的。

参考文献中有s平面映射到z平面的图,方便理解,有需要的同学可以查看。

后向差分映射图

2. 前向差分

前向差分的映射关系为:

从s平面到z平面上来看,其实是将左半平面映射到了z平面上直线 的左边平面上。所以原来在s域上稳定的系统经过前向差分之后的离散系统可能是不稳定的。

前向差分方法映射图

3. Tustin Method(双线性变换方法:Bilinear Transform)

双线性变换方法的映射关系为:

双线性变换方法其实是比较常用的离散化方法啦,毕竟这种方法能够保持映射前后的稳定性是一致的。因为从s平面到z平面的映射来看,该映射将s平面的左半平面映射到z平面的单位圆内。

Tustin方法映射图

4. Tustin Method with pre-warping

双线性变换方法虽然好,但是也有其本身的问题:频偏效应。s域的传递函数的频率和离散化数字传递函数的频率的对应映射关系是非线性的。这种效应可以通过带pre-warping的Tustin方法缓解。频率的对应映射关系为

根据这个频率映射关系,我们就可以针对特定的频率范围附近,也就是我们感兴趣的频率范围附近,将数字滤波器的幅频特性与模拟滤波器的幅频特性相匹配。

5.零极点匹配法

s域的传递函数为:

那么根据零极点匹配得原则,离散化形式如下:

这种方法存在一定的假设前提,同时还有对于离散化得到的数字滤波器加入一个为0的零点来得到一个minimum delay filter这样的技巧,详情可以参见参考。

6. 其他方法

除了上述方法之外呢,还有一些其他方法,例如:

除了上述方法之外呢,还有一些其他方法,例如:

  • 冲激响应不变法
  • 阶跃响应不变法(ZOH法,零阶保持器法)
  • 斜坡响应不变法(一阶保持器法)
  • 幅频响应不变法
  • 相频响应不变法

这些方法都有其特性或有点,在不同的场合下应该有不同的应用。笔者暂时用不到,就暂且不表了。

参考

  1. Vadhavkar P R. MAPPING CONTROLLERS FROM THE S-DOMAIN TO THE Z-DOMAIN USING MAGNITUDE INVARIANCE AND PHASE INVARIANCE METHODS A Thesis by[J]. Wichita State University, 2007.
  2. The design of IIR filters(cont.)
  3. Lecture notes

系统截止频率对应的相角与-180°相差的度数成为相角裕度,为何?

因为所谓相角裕度与增益裕度都是基于频域的控制系统稳定性指标。倘若系统的开环传递函数的截止频率相移为-180°,那么由于负反馈对于正弦波自带180°的相位滞后,那么对于系统而言就相当于正反馈,则系统将处于自激振荡状态。举个简单的例子,前向通道和反馈通道分别为积分环节:

1
2
3
4
fw = tf(1, [1, 0]);
fd = fw ;
sys = feedback(fw, fd);
impulse(sys)

可以看到,系统的脉冲响应是周期为6.28的正弦波,即当存在轻微扰动时,系统就会进入自激振荡状态(因为系统临界稳定);而当系统的输入为0的时候,是输出也为0, 这可以理解为系统无能量输入,所以系统也就不会输出响应。想起模电中运用正反馈构成自激振荡电路的原理,实际工作环境中存在着能量的传播,所以自激振荡电路总能产生振荡的正弦波。理解这一点,对数稳定判据以及奈奎斯特稳定判据就很好理解了。

Impulse response of an unstable system

考虑下图显示的典型的反馈控制系统框图,如果将反馈处切开,那么从给定(Setpoint)到反馈的传递函数实际上就是自控概念中的开环传递函数了。

Controller and Plant

对于正弦信号而言,开环传递函数的幅频特性实际上反映的是反馈回来的信号的放大倍数(也可能是缩小)和相角延迟,这个和系统的稳定性有着密切的关系。如果反馈回来的信号为-180°,并且增益也比较大(大于等于1),那么系统就自然不会稳定。然而在实际系统中,由于功率(能量)的限制,噪声的影响,自激振荡这种状态就很好解释了。

毕竟开环传递函数研究的是反馈相当于系统给定的性质,把握了这句话,电工出身的人所说的环路增益研究系统稳定性就很好理解了。:)