采样信号和傅里叶变换

香农采样定理：采样频率只要大于信号中最高频率的2倍，采样值就可以包含原始信号所有信息

假设如下采样信号

int sampleCount = 10000;
var time = new double[sampleCount];
var xs = new double[sampleCount];
var ys_real = new double[sampleCount];
var ys_img = new double[sampleCount];

//10KHz采样，生成实部数据 和 X轴数据
for(int i=0;i< sampleCount; i++) 
{
	time[i] = i;
	var x= ((double)i) / ((double)sampleCount);
	xs[i] = x;
	ys_real[i] =
		0.5 * Math.Sin(2 * Math.PI * 10 * x ) //叠加10Hz信号
		+ 0.2 * Math.Sin(2 * Math.PI * 80 * x + Math.PI / 4) //叠加80Hz信号
		+ 0.05 * Math.Sin(2 * Math.PI * 200 * x + Math.PI / 6)//叠加200Hz信号
		+ 0.07 * Math.Sin(2 * Math.PI * 400 * x + Math.PI / 8)//叠加400Hz信号
		+ 10;////叠加直流分量
}
//绘制10K个点
_timePlot.Plot.Clear();
_timePlot.Plot.Add.ScatterLine(xs.ToArray(), ys_real.ToArray()); 
_timePlot.Plot.Axes.AutoScale();
_timePlot.Refresh();

//傅里叶向前变换FFT default选型需要导入2的N次方数据 这里选Matlab不需要2的幂数

Fourier.Forward(ys_real, ys_img, MathNet.Numerics.IntegralTransforms.FourierOptions.Matlab);

var result = new double[ys_real.Length];

for (int i = 0; i < result.Count(); i++)
{
	var temp = new Complex32((float)ys_real[i], (float)ys_img[i]);
	result[i] = temp.Magnitude;//复数的模


	//归一化
	//0HZ的直流分量 =振幅/采样点数N 
	if (i == 0) { result[i] = result[i] / sampleCount; }
	else
	{
		//其他==振幅/采样点数N/2
		result[i] = result[i] / (sampleCount / 2);
	}

}
//Scottplot绘制FFT波形
_fftPlot.Plot.Clear();
_fftPlot.Plot.Add.ScatterLine(time.Take(sampleCount/2).ToArray(), result.Take(sampleCount/2).ToArray());
_fftPlot.Plot.Axes.AutoScale();
_fftPlot.Refresh();


频率和整幅都可以匹配

//去除高频部分大于70的频率信号
for (int i = 0; i < ys_real.Length; i++)
{
	if (i > 70)
	{
		ys_real[i] = 0;
		ys_img[i] = 0;
	}
}
//反傅里叶变化
Fourier.Inverse(ys_real, ys_img, FourierOptions.Matlab);
//绘制过滤过的波形 图上黄色波形
_timePlot.Plot.Add.ScatterLine(xs.ToArray(), ys_real.ToArray());
_timePlot.Plot.Axes.AutoScale();
_timePlot.Refresh();

只剩下0.5 * Math.Sin(2 * Math.PI * 10 * x )+ 10;的时域信号

//TIME->FFT-- > FILTER(noise) - IFFT - TIME

用途：

  1. 计算信号的频域信号分析特征
  2. 对信号进行预处理，剔除高频毛刺部分

以上的步骤模拟的频率正好是整数的频次10 80之类 单位是1HZ，现在我们扩展到0.1HZ
采样10K，采样率1K/S


			int N = 10000;//采样10K个
			int fs = 1000;//采样频率
			var n = new double[N];
			var timex = new double[N];
			var ys_real = new double[N];
			var ys_img = new double[N];

			var oneSampleTime = 1d / fs;

			var fftwindow = MathNet.Numerics.Window.Hann(N);
			//生成实部数据 和 X轴数据
			for (int i = 0; i < N; i++)
			{
				n[i] = i;

				timex[i] = oneSampleTime * i;
				ys_real[i] = 
					0//叠加直流分量
					+ 0.5 * Math.Sin(2 * Math.PI * 10.5 * timex[i]) //叠加10Hz信号
					+ 0.2 * Math.Sin(2 * Math.PI * 40.777 * timex[i]) //叠加10Hz信号
					+ 0.1 * Math.Sin(2 * Math.PI * 80 * timex[i] + Math.PI / 4) //叠加80Hz信号
					+ 0.03 * Math.Sin(2 * Math.PI * 200 * timex[i] + Math.PI / 6)//叠加200Hz信号
					+ 0.01 * Math.Sin(2 * Math.PI * 400 * timex[i] + Math.PI / 8)//叠加400Hz信号
					;
			}
			//绘制10K个点
			_timePlot.Plot.Clear();
			_timePlot.Plot.Add.ScatterLine(timex.ToArray(), ys_real.ToArray());
			_timePlot.Plot.Axes.AutoScale();
			_timePlot.Refresh();

			//傅里叶向前变换FFT default选型需要导入2的N次方数据 这里选Matlab不需要2的幂数
			Fourier.Forward(ys_real, ys_img, MathNet.Numerics.IntegralTransforms.FourierOptions.Matlab);

			var result = new double[ys_real.Length];

			for (int i = 0; i < result.Count(); i++)
			{
				var temp = new Complex32((float)ys_real[i], (float)ys_img[i]);
				result[i] = temp.Magnitude;//复数的模


				//归一化
				//0HZ的直流分量 =振幅/采样点数N 
				if (i == 0) { result[i] = result[i] / N; }
				else
				{
					//其他==振幅/采样点数N/2
					result[i] = result[i] / (N / 2);
				}

			}
			var frequenceArray = new double[result.Length];
			for (int i = 0; i < frequenceArray.Length; i++)
			{
				frequenceArray[i] = i * (double)fs / N;
			}
			//Scottplot绘制FFT波形
			_fftPlot.Plot.Clear();
			_fftPlot.Plot.Add.ScatterLine(frequenceArray.ToArray(), result.Take(N / 2).ToArray());
			_fftPlot.Plot.Axes.AutoScale();
			_fftPlot.Refresh();


左图是时域波形
右图是频域波形 注意这次我们插入了两个非整数频率的信号10.5 和 40.777HZ 对40.777的区域放大后可以看到明显的频率泄露，边上的频率也有了振幅
0.5 * Math.Sin(2 * Math.PI * 10.5 * timex[i])
0.2 * Math.Sin(2 * Math.PI * 40.777 * timex[i])

[参考文档]因为离散傅里叶变换DFT只能输出在Fs/N（fs=1000/N=10000）的频率点上的功率，所以当输入频率不在Fs/N的整数倍时，在dft的输出上就没有与输入频率相对应得点(dft输出是离散的)，那么输入频率就会泄漏到所有的输出点上，具体的泄漏分布取决于所采用的窗的连续域复利叶变换
这里带入矩形窗概念，如果我们对生成ys_real[i]数组每个元素都1那就是给时域信号加矩形窗,因为都是1，所以原始的时域信号默认就加了矩形窗
那矩形窗的默认定义就是N的数组，每个元素都是1,ys_real[i]=ys_real[i]fftwindow[i];
现在我对信号采样数N进行扩容1倍2N（N-2N区间的时域值补0） 根据上述原理可以得到0.1/2=0.05HZ的信号

ys_real[i] = 
				0//叠加直流分量
				+ 0.5 * Math.Sin(2 * Math.PI * 10.5 * timex[i]) //叠加10Hz信号
				+ 0.2 * Math.Sin(2 * Math.PI * 40.55 * timex[i]) //叠加10Hz信号
				+ 0.1 * Math.Sin(2 * Math.PI * 80 * timex[i] + Math.PI / 4) //叠加80Hz信号
				+ 0.03 * Math.Sin(2 * Math.PI * 200.5 * timex[i] + Math.PI / 6)//叠加200Hz信号
				+ 0.01 * Math.Sin(2 * Math.PI * 400 * timex[i] + Math.PI / 8)//叠加400Hz信号
				;
//补零
var zero = new double[N * 2];
Array.Copy(ys_real,zero, ys_real.Length);
ys_real = zero;
ys_img = new double[N * 2];


观测到40.55Hz的频率可以抓到但是两边都有泄露，那我们需要加个hann窗什么的观察一下改善

int N = 10000;//采样10K个
int fs = 1000;//采样频率
var n = new double[N];
var timex = new double[N];
var ys_real = new double[N];
var ys_img = new double[N];

var oneSampleTime = 1d / fs;

var fftwindow = MathNet.Numerics.Window.Hann(N);
//生成实部数据 和 X轴数据
for (int i = 0; i < N; i++)
{
	n[i] = i;

	timex[i] = oneSampleTime * i;
	ys_real[i] = 
		0//叠加直流分量
		+ 0.5 * Math.Sin(2 * Math.PI * 10.5 * timex[i]) //叠加10Hz信号
		+ 0.2 * Math.Sin(2 * Math.PI * 40.55 * timex[i]) //叠加10Hz信号
		+ 0.1 * Math.Sin(2 * Math.PI * 80 * timex[i] + Math.PI / 4) //叠加80Hz信号
		+ 0.03 * Math.Sin(2 * Math.PI * 200.5 * timex[i] + Math.PI / 6)//叠加200Hz信号
		+ 0.01 * Math.Sin(2 * Math.PI * 400 * timex[i] + Math.PI / 8)//叠加400Hz信号
		;
	//先加窗后补0
	if (i < N)
	{
		ys_real[i] = ys_real[i] * fftwindow[i ]; //时域信号的乘积=频域信号的卷积
	}
	//补零
	var zero = new double[N * 2];
	Array.Copy(ys_real,zero, ys_real.Length);
	ys_real = zero;
	ys_img = new double[N * 2];
	//if (i > N / 2) { ys_real[i] = 0; }

}
//绘制10K个点
_timePlot.Plot.Clear();
_timePlot.Plot.Add.ScatterLine(timex.ToArray(), ys_real.ToArray());
_timePlot.Plot.Axes.AutoScale();
_timePlot.Refresh();

//傅里叶向前变换FFT default选型需要导入2的N次方数据 这里选Matlab不需要2的幂数
Fourier.Forward(ys_real, ys_img, MathNet.Numerics.IntegralTransforms.FourierOptions.Matlab);

var result = new double[ys_real.Length];

for (int i = 0; i < result.Count(); i++)
{
	var temp = new Complex32((float)ys_real[i], (float)ys_img[i]);
	result[i] = temp.Magnitude;//复数的模


	//归一化
	//0HZ的直流分量 =振幅/采样点数N 
	if (i == 0) { result[i] = result[i] / N; }
	else
	{
		//其他==振幅/采样点数N/2
		result[i] = result[i] / (N / 2);
	}

}
var frequenceArray = new double[result.Length];
for (int i = 0; i < frequenceArray.Length; i++)
{
	frequenceArray[i] = i * (double)fs / (2*N);
}
//Scottplot绘制FFT波形
_fftPlot.Plot.Clear();
_fftPlot.Plot.Add.ScatterLine(frequenceArray.ToArray(), result.Take(N / 2).ToArray());
_fftPlot.Plot.Axes.AutoScale();
_fftPlot.Refresh();


加窗后的时域波形N-2N区间的没有绘制出来就是一根直线；

可以观察到泄露小了很多；
思路就是采集一段时间长度的信号，如果要观察的频率有小数点几位就要适当的进行扩容（补0），扩容前加窗函数抑制泄露

参考文档
[傅里叶变化-Math.Net中Fourier变换使用详解]
https://blog.csdn.net/u011555996/article/details/136403452
[如何 FFT(快速傅里叶变换) 求幅度、频率（超详细 含推导过程）]
https://blog.csdn.net/weixin_39591031/article/details/110392352
[信号处理之快速傅里叶变换(FFT)-通俗易懂]
https://blog.csdn.net/zhang0558/article/details/136840352

posted @   YYAN1987  阅读(66)  评论(0)  编辑  收藏  举报
相关博文：
· WPFLocalizeExtension 多语言组件使用和原理及拓展
· Matlab 实现抽样定理
· 奈奎斯特采样定理
· matlab FFT -2020
· 傅里叶变换
阅读排行：
· 被坑几百块钱后，我竟然真的恢复了删除的微信聊天记录！
· 没有Manus邀请码？试试免邀请码的MGX或者开源的OpenManus吧
· 【自荐】一款简洁、开源的在线白板工具 Drawnix
· 园子的第一款AI主题卫衣上架——"HELLO! HOW CAN I ASSIST YOU TODAY
· 无需6万激活码！GitHub神秘组织3小时极速复刻Manus，手把手教你使用OpenManus搭建本
点击右上角即可分享
微信分享提示