我们为什么用卷积? Why should we use convolution?

Why should we use convolution?



问题限定: 仅对离散信号做分析


首先要回答什么是卷积的问题.


              In mathematics and, in particular, functional analysisconvolution is a mathematical operation on two functions f and g, producing a third function that is typically viewed as a modified version of one of the original functions, giving the area overlap between the two functions as a function of the amount that one of the original functions is translated


                                                                                                                          —— wikipedia


对于卷积这个东东,一般看到这两个字的人就头大,往往望而生畏.

上述定义严谨,而《The sciencetist and engineer's guide to DSP》有更为简单直接的描述


Convolution is a mathematical way of combining two signals to form a third signal.


                 Convolution is a formal mathematical operation, just as multiplication, addition, and integration. Addition takes two numbers and produces a third number, while convolution takes two signals and produces a third signal.

不要有太大的心理压力,它不过就是一种数学计算操作而已,用两个输入信号合成第三个信号输出.至于具体的怎么操作,后面一步步讲清楚.


先抛出个问题,然后借助解答这个问题去回答“what is convolution ?".


这里有两组很简单信号x[n], h[n],这两组信号做卷积处理之后得到y[n]


问题来了,x[n] 和 h[n]怎么操作就可以得到y[n]呢?


一步步分解

第一步,图1,x[0]h[n-0] ,这里表达的就是把x[0]依次 和h[0]~h[3]做乘法,把结果输出到此时对应的output[n](图中的方块)

第二步,图2,x[1]h[n-1] ,这里表达的就是把x[1]依次 和h[0]~h[3]做乘法,把结果输出到此时对应的output[n](图中的方块)

第三步,图3,x[2]h[n-2] ,这里表达的就是把x[2]依次 和h[0]~h[3]做乘法,把结果输出到此时对应的output[n](图中的方块)

... ...

第九步,图9,x[8]h[n-8] ,这里表达的就是把x[8]依次 和h[0]~h[3]做乘法,把结果输出到此时对应的output[n](图中的方块)


So ? 看出规律了?

之前我故意没有给出卷积的公式定义

现在再看卷积的公式定义一点清楚得多!


用Octave(或者 matlab)实现一个demo

利用了这幅图的数据



%code writer	:	EOF
%code date	:	2014.09.11
%e-mail		:	jasonleaster@gmail.com
%code file	:	convolution_input_side.m
%code purpose:
%             A demo for convolution.
clear all
close all

x = [0 -1 -1.2 2 1.4 1.4 0.8 0 -0.8];
h = [1 -0.6 -0.4 -0.1];

length_x = size(x,2);
length_h = size(h,2);

output = zeros(1,length_h+length_x);

figure(1);
subplot(121);
scatter(1:length_x,x,'r','filled');
title('x[n]');
subplot(122);
scatter(1:length_h,h,'g','filled');
title('h[n]');

for tmp_1 = 1:length_x

	for tmp_2 = 1:length_h
	
		output(tmp_1+tmp_2-1) = output(tmp_1+tmp_2-1) + x(tmp_1)*h(tmp_2);
	end
end

figure(2);
scatter(1:size(output,2),output,'filled');
title('output[n]');



到此,回答了什么是卷积.或者说更深刻体会了什么是卷积


为什么我们要用卷积这种东东?


 

Why should we use convolution? God said I just want to make filter more simple. : )


时域中的卷积就对应于频域中的乘积,那么频域中的卷积对应于时域中的乘积


卷积为滤波而生.


这里是两个利用卷积做滤波器的例子



图a是低通滤波,图b是高通滤波,一开始我也纳闷Input signal确实复合了多个频率的信号,怎么耶和中间的信号做卷积之后就有滤波的效果了呢?


低通滤波做卷积的这个信号有什么特点吗??

filter = [0 0.1 0.16 0.25 0.25 0.16 0.1 0]*0.125;
scatter(1:8,real(fft(filter)),'filled');

假设图中的数据是filter(我自己模拟的数据)

可以用频谱分析发现

这货的频谱是低频部分高的!抑制高频部分


把两个输入信号都变换到频率领域里面去,然后想乘,再把结果逆变换到时间领域不久得到低通滤波的结果了么!



这里再给出一种低通滤波的demo

首先,利用一下两个信号合成一个输入信号,模拟带高频噪声的输入


利用上面两个信号,我模拟合成了下面第一幅图的信号,作为带高频噪声的输入信号,而下面的第二幅图就是我们的做卷积的另外一个信号(低通滤波就靠它了!)


将以上两个信号做卷积得到输出信号

漂亮!几乎完美的滤去的高频噪声.



实现:

%code writer	:	EOF
%code date	:	2014.09.11
%e-mail		  :	jasonleaster@gmail.com
%code file	:	convolution_input_side.m
%code purpose:
%        Just a demo for low-pass filter by convolution.


clear all
close all

% 'y' is a simulation for original pure signal
y = sin(0.5*pi*[0:0.05:8]); 
% 'u' is a simulation for noise.
u = 0.2*sin(8*pi*[0:0.05:8]);

figure(1);
subplot(211);
plot(1:size(y,2),y);
title('original pure signal'); %BUG, I don't know why I can't write title on this figure.
subplot(212);
plot(1:size(u,2),u);
title('noise');

input_signal = y + u;

x = input_signal;
h = [1 1 1 1 1 1 1 1 1 1]*0.1;

length_x = size(x,2);
length_h = size(h,2);

output_signal = zeros(1,length_h+length_x);

figure(2);
subplot(211);
plot(1:length_x,x);
subplot(212);
scatter(1:length_h,h,'filled');

for tmp_1 = 1:length_x

	for tmp_2 = 1:length_h
	
		output_signal(tmp_1+tmp_2-1) = output_signal(tmp_1+tmp_2-1) + x(tmp_1)*h(tmp_2);
	end
end

figure(3);
plot(1:size(output_signal,2),output_signal);
title('Output signal');




哈哈是不是想看看这个家伙的频谱图呢?



怎么它就可以低通滤波捏?下面给出了上面那个家伙的频谱图,看!低通滤波!只有在低频地带有赋值为1,其他的都是0!



当我们想要什么类型的滤波器,我们可以先构想频谱表达,然后通过逆向傅立叶变换(实际上是逆向FFT),然后得到和输入信号做卷积的因子,然后做卷积,即可实现滤波!



 

Why should we use convolution? God said I just want to make filter more simple. : )


update: 2014.11.13

给出卷积的C语言实现demo

/*********************************************************
code writer 	: EOF
code file   	: conv_demo.c
code date   	: 2014.11.13
e-mail		: jasonleaster@gmail.com

code description:

	A implementation of convolution in C lanuage

**********************************************************/
#include <stdio.h>
#include <stdlib.h>

#define DEBUG

int* init_input(int number)
{
	if(number <= 0)
	{
		printf("number < 0! in function %s\n",__FUNCTION__);
		return NULL;
	}

	int* p_input = (int*)malloc(sizeof(int)*number);

	if(!p_input)
	{
		printf("malloc error!\n");

		return NULL;
	}

	int temp = 0;
	for(temp = 0;temp < number;temp++)
	{
		if(scanf("%d",(p_input+temp)) < 0)
		{
			return NULL;
		}
	}

	return p_input;
}

int* init_unit_responce(int number)
{
	if(number <= 0)
	{
		printf("number < 0! in function %s\n",__FUNCTION__);
		return NULL;
	}

	int* p_unit_responce = (int*)malloc(sizeof(int)*number);

	if(!p_unit_responce)
	{
		printf("malloc error!\n");
		return NULL;
	}

	int temp = 0;
	for(temp = 0;temp < number;temp++)
	{
		if(scanf("%d",(p_unit_responce +temp)) < 0)
		{
			return NULL;
		}
	}

	return p_unit_responce;
}

void release(int* p_input,int* p_unit_responce)
{
	free(p_input);
	free(p_unit_responce);
}

int* conv(int* p_input,int ip_num,int* p_unit_responce,int ur_num)
{
	if(!p_input || !p_unit_responce || ip_num <= 0 || ur_num <= 0)
	{
		printf("Data losted!\n");

		return NULL;
	}

	int* p_output = (int*)malloc(sizeof(int)*(ip_num + ur_num - 1));

	if(!p_output)
	{
		printf("malloc error! in function %s\n",__FUNCTION__);
		return NULL;
	}

	int tmp = 0;
	int idx = 0;
	int sum = 0;

	for(tmp = 0;tmp < ip_num + ur_num -1;tmp++)
	{
		sum = 0;
		for(idx = 0; (tmp - idx) >= 0 && idx < ur_num;idx++)
		{
			if((tmp -idx) < ip_num)
			{
				sum += p_unit_responce[idx]*p_input[tmp-idx];
			}
		}

		*(p_output + tmp) = sum;

		#ifdef DEBUG
			printf("\t%d",sum);
		#endif
	}

	#ifdef DEBUG
		printf("\n");
	#endif

	return p_output;
}

int main()
{
	int ip_num = 0;

	printf("Hey guys! Just input the size of input-data!\n");

	while(!scanf("%d",&ip_num))
	{
		printf("input data type error! Input again!\n");

		getchar();// clear the '\n' after we inputed.
	}

	int* p_input = NULL;	

	p_input = init_input(ip_num);

	if(!p_input)
	{
		goto out;
	}

	printf("What about the size of unit-responce!\n");

	int ur_num = 0;
	while(!scanf("%d",&ur_num))
	{
		printf("input data type error! Input again!\n");

		getchar();// clear the '\n' after we inputed.
	}

	int* p_unit_responce = NULL;

	p_unit_responce = init_unit_responce(ur_num);

	if(!p_unit_responce)
	{
		goto out;
	}

	conv(p_input,ip_num,p_unit_responce,ur_num);

out:
	release(p_input,p_unit_responce);
	return 0;

}

测试结果:


用matlab的conv函数来检验我们程序的正确性,对比可知,我们程序实现是正确的。



评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值