傅里叶变换的例题

千次阅读

2022-02-12 23:35:10

目的：去除噪音。时间和信号的关系是：gaussNoisy(t) （g(t)）

noiseAmp：控制噪音的大小。noiseAmp为0时图像会变成光滑的信号。

import numpy as np

def gaussNoisy(t, noiseAmp):
    noise = noiseAmp*(np.random.randn(len(t)))
    return np.exp(-0.5*t*t) * (1.0+noise)

N = 256
t = np.linspace(-4.0, 4.0, N)
g = gaussNoisy(t, 0.1)

（閾値は試行錯誤により決定しなさい：范围根据错误决定，我试了好多次xlim和ylim找噪音）

以下为解决方法的程序：

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as scifft
plt.close('all')
def gaussNoisy(t, noiseAmp):
    noise = noiseAmp*(np.random.randn(len(t)))
    return np.exp(-0.5*t*t) * (1.0+noise)

N = 256
t = np.linspace(-4.0, 4.0, N)
g = gaussNoisy(t, 0.1)


fig1, ax1 = plt.subplots(2, 1, figsize=(9, 7))
ax1[0].plot(t,g, color='red') 
ax1[0].set_xlim(-4, 4) 
ax1[0].set_ylim(0,1.3) 
ax1[0].set_xlabel('t') 
ax1[0].set_ylabel('g(t)')


num_data =len(t)
freq = scifft.fftfreq(num_data, d=0.03125)  
freq = scifft.fftshift(freq)

f_hat = scifft.fft(g,n=num_data)
f_hat = scifft.fftshift(f_hat) 
f_spect = f_hat * np.conj(f_hat)/num_data
fig1, ax2 = plt.subplots(4, 1, figsize=(9, 7))
ax2[0].plot(freq, np.real(f_hat), color='blue')
ax2[0].set_xlim(-128,128)
ax2[0].set_ylim(-1.2, 1.2) 
ax2[0].set_xlabel('freq (Hz)') 
ax2[0].set_ylabel('Re[f_hat]') 

# f_hat(f)作图（虚部imag）
ax2[1].plot(freq, np.imag(f_hat), color='blue')
ax2[1].set_xlim(-128,128) 
ax2[1].set_ylim(-1.2, 1.2) 
ax2[1].set_xlabel('freq (Hz)') 
ax2[1].set_ylabel('Im[f_hat]')
 
# f_hat的spect作图
ax2[2].plot(freq,f_spect, color='blue',label="before filtering") #傅里叶变换：可以看出有噪音
ax2[2].set_xlim(-128,128) 
ax2[2].set_ylim(-0.01, 0.01) 
ax2[2].set_xlabel('freq (Hz)') 
ax2[2].set_ylabel('f_spect') 
ax2[2].legend()

indices = f_spect>0.05 
#print("indices=", indices)
 

f_spect_filtered = f_spect * indices
f_hat_filtered = f_hat * indices
 
ax2[3].plot(freq,f_spect_filtered, color='red',label="after filtering")
ax2[3].set_xlim(-128,128)
ax2[3].set_ylim(-0.05, 0.05) 
ax2[3].set_xlabel('freq (Hz)') 
ax2[3].set_ylabel('f_spect(freq)') 
ax2[3].legend()
 
#逆傅里叶变换准备：要用fftshift
f_hat_filtered = scifft.fftshift(f_hat_filtered)
f_filterd = scifft.ifft(f_hat_filtered,n=num_data)
 
ax1[1].set_xlim(-4, 4) 
ax1[1].set_ylim(0,1.3) 
ax1[1].plot(t, np.real(f_filterd), color='blue') 
ax1[1].set_xlabel('time (sec)') 
ax1[1].set_ylabel('f_filtered(t)') 
 
plt.show()
plt.tight_layout()

python

更多相关内容

傅里叶变换例题
2014-05-28 19:10:40

常用的题型，很有用。主要飚速了傅里叶变换，原理作用等等

收起
信号与系统实验傅里叶变换.docx
2022-03-22 16:15:05

信号与系统实验傅里叶变换

收起
课时17 傅里叶变换法求解微分方程.pdf
2021-12-13 21:38:43

数学物理方程Mathematical Physics任永志物理与光电工程学院哈尔滨工程大学renyongzhi@hrbeu.edu.cn1( )( )2ikxf kf x edx1( )( )2ikxf xf kdxe 傅里叶变换定义傅里叶变换复习常用...

收起

二维离散傅里叶变换——数字图像处理学习六（C++版）

千次阅读 2021-07-06 18:12:26

参考资料视频：https://www.bilibili.com/video/BV1za411F76U?from=search&seid=15738765176302809882 二、二维离散傅里叶变换的计算过程例题：解：根据二维离散傅里叶计算公式其中M=4,N=4 和欧拉公式可得： ...

一、名词解释

FT(Fourier Transform),傅里叶变换

傅立叶变换，表示能将满足一定条件的某个函数表示成三角函数或者它们的积分的线性组合。在不同的研究领域，傅立叶变换具有多种不同的变体形式，如连续傅立叶变换和离散傅立叶变换。

1)傅里叶变换

$\large F\left ( \omega \right )= F\left [ f\left ( t \right ) \right ]= \int_{-\infty }^{+\infty }f\left ( t \right )e^{-i\omega t}dt$

2)傅里叶逆变换

$\large f\left (t \right )= F^{-1}\left [ F\left ( \omega \right ) \right ]=\frac{1}{2\pi } \int_{-\infty }^{+\infty }F\left (\omega \right )e^{i\omega t}d\omega$

DFT(Discrete Fourier Transform),离散傅里叶变换

1)二维离散傅里叶变换

$\large F\left ( u,v \right )=\sum_{x= 0}^{M-1}\sum_{y= 0}^{N-1}f\left ( x,y \right )e^{-j2\pi \left ( \frac{ux}{M} +\frac{vy}{N}\right )}$

2)二维离散傅里叶拟变换

$\large f\left ( x,y \right )=\frac{1}{MN}\sum_{u= 0}^{M-1}\sum_{v= 0}^{N-1}F\left ( u,v \right )e^{j2\pi \left ( \frac{ux}{M} +\frac{vy}{N}\right )}$

FFT(fast Fourier transform),快速傅里叶变换

快速傅里叶变换是计算机计算离散傅里叶变换的高效、快速的计算方法。采用这种算法能使计算机计算离散傅里叶变换需要的时间复杂度减少。

参考资料视频：https://www.bilibili.com/video/BV1za411F76U?from=search&seid=15738765176302809882

二、二维离散傅里叶变换的计算过程

例题： $e^{jx}= cosx+jsinx$

$f\left ( x,y \right )=\begin{bmatrix} 1 & 7& 8& 9\\ 5& 2& 6& 7\\ 3& 4& 7& 8\\ 2& 1& 4& 7 \end{bmatrix}$

解：根据二维离散傅里叶计算公式

$\large F\left ( u,v \right )=\sum_{x= 0}^{M-1}\sum_{y= 0}^{N-1}f\left ( x,y \right )e^{-j2\pi \left ( \frac{ux}{M} +\frac{vy}{N}\right )}$

其中M=4,N=4

和欧拉公式

$\large e^{jx}= cosx+jsinx$

可得：

其中 $\large u=1,2,3,4 v=1,2,3,4 x=1,2,3,4 y=1,2,3,4$

则当 $\large u=0,v=0$ 时

$\large F\left ( 0,0 \right )=f\left ( 0,0 \right )+f\left ( 0,1 \right )+\cdots +f\left ( 3,3 \right )=1+7+8+9+\cdots +4+7=81$

为了加快计算，利用代码进行计算

#define PI 3.1415926
int main()
{
	long double f[16] = { 1.0,7.0,8.0,9.0,5.0,2.0,6.0,7.0,3.0,4.0,7.0,8.0,2.0,1.0,4.0,7.0 };	
long double dft2_data_real[16] = { 0 };
	long double dft2_data_xu[16] = { 0 };
	long double fixed_factor_for_axisX = (-2 * PI) / 4;//-2π/M
	long double fixed_factor_for_axisY = (-2 * PI) / 4;//-2π/N
	for (int u = 0; u < 4; u++)
	{
		for (int v = 0; v < 4; v++)
		{
			for (int x = 0; x < 4; x++)
			{
				for (int y = 0; y < 4; y++)
				{
					long double powerX = u * x * fixed_factor_for_axisX; //-2πux/M
					long double powerY = v * y * fixed_factor_for_axisY; //-2πvy/N
					long double real = f[y + x * 4] * cos(powerX + powerY);
					long double xu = f[y + x * 4] * sin(powerX + powerY);
					dft2_data_real[u * 4 + v] = dft2_data_real[u * 4 + v]+real;
					dft2_data_xu[u * 4 + v] = dft2_data_xu[u * 4 + v]+xu;
				}
			}
		}
	}
	for (int u = 0; u < 4; u++)
	{
		for (int v = 0; v < 4; v++)
		{
			cout << "real:" << dft2_data_real[u * 4 + v] << "xu:" << dft2_data_xu[u * 4 + v] << endl;
		}
	}
}

执行结果

所以结果为

$\large F\left ( u,v \right )=\begin{bmatrix} 81& -14+17j& -9& -14-17j\\ 3-6j& -4-3j& -5-4j& -2+j\\ 13& -8-5j& -9& -8+5j\\ 3+6j& -2& -5+4j& -4+3j\end{bmatrix}$

三、图像的离散傅里叶变换

原灰度图：

展示频谱图

void dft()
{
CPaintDC dc(this);

	FILE* fp1;
	if (fopen_s(&fp1, "C:\\Users\\Administrator\\Desktop\\BMP1\\BMP1\\1.bmp", "rb") != 0) {
		return;
	}
	bmpFileHeader bfHead1;
	bmpInfoHeader biHead1;
	fread_s(&bfHead1, 14, 14, 1, fp1);
	fread_s(&biHead1, 40, 40, 1, fp1);

	uchar* p1 = (uchar*)malloc(biHead1.biSizeImage * sizeof(uchar));
	fseek(fp1, bfHead1.bfOffBits, SEEK_SET);
	fread_s(p1, biHead1.biSizeImage * sizeof(uchar), biHead1.biSizeImage * sizeof(uchar), 1, fp1);
	///灰度化
	unsigned long x = 0;
	uchar* gray = (uchar*)malloc(biHead1.biHeight * biHead1.biWidth * sizeof(uchar));
	for (int i = 0; i < biHead1.biHeight; i++)
	{
		for (int j = 0; j < biHead1.biWidth; j++)
		{
			gray[i * biHead1.biWidth + j] = 0.2989 * p1[x + 2] + 0.5870 * p1[x + 1] + 0.1140 * p1[x];
			x += 3;
		}
	}
	//转成double型
	double* double_data = new double[biHead1.biHeight * biHead1.biWidth];
	for (int i = 0; i < biHead1.biHeight; i++)
		for (int j = 0; j < biHead1.biWidth; j++)
		{
			double_data[i * biHead1.biWidth + j] = (double)gray[i * biHead1.biWidth + j];
		}
	// 对 double_data 进行傅里叶变换
	double* dft2_data_real = new double[biHead1.biHeight * biHead1.biWidth]();//初始化0
	double* dft2_data_xu = new double[biHead1.biHeight * biHead1.biWidth]();//初始化0
	double fixed_factor_for_axisX = (-2 * PI) / biHead1.biHeight;//-2π/M
	double fixed_factor_for_axisY = (-2 * PI) / biHead1.biWidth;//-2π/N
	for (int u = 0; u < biHead1.biHeight; u++)
	{
		for (int v = 0; v < biHead1.biWidth; v++)
		{
			for (int x = 0; x < biHead1.biHeight; x++)
			{
				for (int y = 0; y < biHead1.biWidth; y++)
				{
					double powerX = u * x * fixed_factor_for_axisX; //-2πux/M
					double powerY = v * y * fixed_factor_for_axisY; //-2πvy/N
					double real = double_data[y + x * biHead1.biWidth] * cos(powerX + powerY);
					double xu = double_data[y + x * biHead1.biWidth] * sin(powerX + powerY);
					dft2_data_real[u * biHead1.biWidth + v] = dft2_data_real[u * biHead1.biWidth + v] + real;
					dft2_data_xu[u * biHead1.biWidth + v] = dft2_data_xu[u * biHead1.biWidth + v] + xu;
				}
			}
		}
	}

	// 将傅里叶数组归一化
	// 取模
	double* mold = new double[biHead1.biHeight * biHead1.biWidth];
	for (int i = 0; i < biHead1.biHeight * biHead1.biWidth; i++)
	{
		mold[i] = sqrt(dft2_data_real[i] * dft2_data_real[i] + dft2_data_xu[i] * dft2_data_xu[i]);
	}
	// 获取最小值
	double min = mold[0];
	for (int i = 0; i < biHead1.biHeight * biHead1.biWidth; i++)
	{
		if (mold[i] < min)
			min = mold[i];
	}
	// 获取去掉前几大值的最大值
	double maxqueue[20] = { 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0. }, max;
	for (int i = 0; i < biHead1.biHeight * biHead1.biWidth; i++) 
	{
		if (mold[i] > maxqueue[0])
			maxqueue[0] = mold[i];
	}
	for (int j = 1; j < 20; j++) 
	{
		for (int i = 0; i < biHead1.biHeight * biHead1.biWidth; i++) 
		{
			if (mold[i] > maxqueue[j] && mold[i] < maxqueue[j - 1])
				maxqueue[j] = mold[i];
		}
	}
	max = maxqueue[19];
	unsigned char* normalized_data = new unsigned char[biHead1.biHeight * biHead1.biWidth];
	for (int i = 0; i < biHead1.biHeight; i++)
		for (int j = 0; j < biHead1.biWidth; j++)
		{
			unsigned char t = (unsigned char)((mold[i * biHead1.biWidth + j] - min) / (max - min) * 255);
			if (t > 255)
				t = 255;
			normalized_data[i * biHead1.biWidth + j] = t;
		}

	// 将图像中心化
	unsigned char* center_data = new unsigned char[biHead1.biHeight * biHead1.biWidth];
	for (int u = 0; u < biHead1.biHeight; u++)
	{
		for (int v = 0; v < biHead1.biWidth; v++) 
		{
			if ((u < (biHead1.biHeight/ 2)) && (v < (biHead1.biWidth / 2))) 
				center_data[v + u * biHead1.biWidth] =normalized_data[biHead1.biWidth / 2 + v + (biHead1.biHeight / 2 + u) * biHead1.biWidth];
			else if ((u < (biHead1.biHeight / 2)) && (v >= (biHead1.biWidth / 2)))
				center_data[v + u * biHead1.biWidth] =normalized_data[(v - biHead1.biWidth / 2) + (biHead1.biHeight / 2 + u) * biHead1.biWidth];
			else if ((u >= (biHead1.biHeight / 2)) && (v < (biHead1.biWidth / 2)))
				center_data[v + u * biHead1.biWidth] =normalized_data[(biHead1.biWidth / 2 + v) + (u - biHead1.biHeight / 2) * biHead1.biWidth];
			else if ((u >= (biHead1.biHeight / 2)) && (v >= (biHead1.biWidth / 2)))
				center_data[v + u * biHead1.biWidth] =normalized_data[(v - biHead1.biWidth / 2) + (u - biHead1.biHeight / 2) * biHead1.biWidth];
		}
	}
	// 向中心化的数组填充空间，频谱图显示
	int bytesPerLine = (biHead1.biWidth * 8 + 31) / 32 * 4;
	unsigned char* dst_data = new unsigned char[bytesPerLine * biHead1.biHeight];
	for (int i = 0; i < biHead1.biHeight; i++)
		for (int j = 0; j < biHead1.biWidth; j++)
		{
			dst_data[i * bytesPerLine + j] = center_data[i * biHead1.biWidth + j];
			dc.SetPixel(j/time, (biHead1.biHeight - i)/time, RGB(dst_data[j + i * biHead1.biWidth], dst_data[j + i * biHead1.biWidth], dst_data[j + i * biHead1.biWidth]));
		}

}

频谱图显示结果：

离散傅里叶逆变换：

void ni()
{
CPaintDC dc(this);

	FILE* fp1;
	if (fopen_s(&fp1, "C:\\Users\\Administrator\\Desktop\\BMP1\\BMP1\\1.bmp", "rb") != 0) {
		return;
	}
	bmpFileHeader bfHead1;
	bmpInfoHeader biHead1;
	fread_s(&bfHead1, 14, 14, 1, fp1);
	fread_s(&biHead1, 40, 40, 1, fp1);

	uchar* p1 = (uchar*)malloc(biHead1.biSizeImage * sizeof(uchar));
	fseek(fp1, bfHead1.bfOffBits, SEEK_SET);
	fread_s(p1, biHead1.biSizeImage * sizeof(uchar), biHead1.biSizeImage * sizeof(uchar), 1, fp1);
	///灰度化
	unsigned long x = 0;
	uchar* gray = (uchar*)malloc(biHead1.biHeight * biHead1.biWidth * sizeof(uchar));
	for (int i = 0; i < biHead1.biHeight; i++)
	{
		for (int j = 0; j < biHead1.biWidth; j++)
		{
			gray[i * biHead1.biWidth + j] = 0.2989 * p1[x + 2] + 0.5870 * p1[x + 1] + 0.1140 * p1[x];
			x += 3;
		}
	}
	//转成double型
	double* double_data = new double[biHead1.biHeight * biHead1.biWidth];
	for (int i = 0; i < biHead1.biHeight; i++)
		for (int j = 0; j < biHead1.biWidth; j++)
		{
			double_data[i * biHead1.biWidth + j] = (double)gray[i * biHead1.biWidth + j];
		}
	// 对 double_data 进行傅里叶变换
	double* dft2_data_real = new double[biHead1.biHeight * biHead1.biWidth]();//初始化0
	double* dft2_data_xu = new double[biHead1.biHeight * biHead1.biWidth]();//初始化0
	double fixed_factor_for_axisX = (-2 * PI) / biHead1.biHeight;//-2π/M
	double fixed_factor_for_axisY = (-2 * PI) / biHead1.biWidth;//-2π/N
	for (int u = 0; u < biHead1.biHeight; u++)
	{
		for (int v = 0; v < biHead1.biWidth; v++)
		{
			for (int x = 0; x < biHead1.biHeight; x++)
			{
				for (int y = 0; y < biHead1.biWidth; y++)
				{
					double powerX = u * x * fixed_factor_for_axisX; //-2πux/M
					double powerY = v * y * fixed_factor_for_axisY; //-2πvy/N
					double real = double_data[y + x * biHead1.biWidth] * cos(powerX + powerY);
					double xu = double_data[y + x * biHead1.biWidth] * sin(powerX + powerY);
					dft2_data_real[u * biHead1.biWidth + v] = dft2_data_real[u * biHead1.biWidth + v] + real;
					dft2_data_xu[u * biHead1.biWidth + v] = dft2_data_xu[u * biHead1.biWidth + v] + xu;
				}
			}
		}
	}
	//逆傅里叶变换
//double* double_data = new double[biHead1.biHeight * biHead1.biWidth];
	fixed_factor_for_axisX = (2 * PI) / biHead1.biHeight;//2π/M
	fixed_factor_for_axisY = (2 * PI) / biHead1.biWidth;//2π/N
	double_data[0] = 0;
	for (int x = 0; x < biHead1.biHeight; x++)
	{
		for (int y = 0; y < biHead1.biWidth; y++)
		{
			for (int u = 0; u < biHead1.biHeight; u++)
			{
				for (int v = 0; v < biHead1.biWidth; v++)
				{
					double powerX = u * x * fixed_factor_for_axisX; //2πux/M
					double powerY = v * y * fixed_factor_for_axisY; //2πvy/N
					double real = dft2_data_real[v + u * biHead1.biWidth] * cos(powerX + powerY) - dft2_data_xu[v + u * biHead1.biWidth] * sin(powerX + powerY);
					double_data[x * biHead1.biWidth + y] = double_data[x * biHead1.biWidth + y] + real;
				}
			}
			double_data[x * biHead1.biWidth + y] = double_data[x * biHead1.biWidth + y] / (biHead1.biHeight * biHead1.biWidth);
		}
	}
//逆傅里叶变换图片显示
	for (int i = 0; i < biHead1.biHeight; i++)
		for (int j = 0; j < biHead1.biWidth; j++)
		{
			dc.SetPixel(j / time, (biHead1.biHeight - i) / time, RGB(double_data[j + i * biHead1.biWidth], double_data[j + i * biHead1.biWidth], double_data[j + i * biHead1.biWidth]));
		}


}

逆变换结果：

注意：本文处理图像时未利用FFT方法，在实际应用中，处理M*N像素的图片DFT方法时间复杂度为M*N*M*N，FFT作者还未清楚，待日后学习补充

收起

展开全文

图像处理

matlab中离散化傅里叶变换,Matlab模拟傅里叶变换
2021-04-21 19:31:17

傅里叶变换是我们最早开始接触的时频域变换方法，虽然经常使用，知道怎么用纸笔计算，但是还从来没有在电脑中模拟过，正好现在开始学习数字信号处理，借着这个机会再学习如何在电脑上模拟傅里叶变换。以下大部分内容...

傅里叶变换是我们最早开始接触的时频域变换方法，虽然经常使用，知道怎么用纸笔计算，但是还从来没有在电脑中模拟过，正好现在开始学习数字信号处理，借着这个机会再学习如何在电脑上模拟傅里叶变换。

以下大部分内容来自Digital Signal Processing Using Matlab和数字信号处理教程程佩青

此次选择的软件平台为Matlab。

由于Matlab无法处理无限长序列，所以需要处理的信号必须是有限长的。

连续时间傅里叶变换

傅里叶变换的公式为：

为了在计算机中模拟傅里叶变换，我们将积分变为求和的方式，上下限也从正无穷到负无穷变为一段长度M，dt需要尽可能小

在Matlab中，函数的自变量因变量的集合都是使用矩阵来存储的，从矩阵的角度来看傅里叶变换的公式如下：

角频率向量定义为

时间向量定义为

因此矩阵指数可写为

整个傅里叶变换可写为

Xa = xa * exp(-1j*t'*W) * Dt;

具体实现

其实下面这个例子是Digital Signal Processing Using Matlab中的，来自P64页，不过想到都看到这里了还要读者翻书不太好，就一起放上来了。

定义

先进行数学上的分析，

MATLAB实现如下：

% Analog Signal

Dt = 0.00005;

t = -0.005:Dt:0.005;

xa = exp(-1000*abs(t));

% Continuous-time Fourier Transform

Wmax = 2*pi*2000;

K = 500;

k = 0:1:K;

W = k*Wmax/K;

Xa = xa * exp(-1j*t'*W) * Dt;

Xa = abs(Xa);

W = [-fliplr(W), W(2:501)];

Xa = [fliplr(Xa), Xa(2:501)];

subplot(2,1,1);

plot(t*1000,xa);

xlabel('t in msec.');

ylabel('xa(t)');

title('Analog Signal');

subplot(2,1,2);

plot(W/(2*pi*1000),Xa*1000);

xlabel('Frequency in KHz'); ylabel('Xa(jW)*1000');

title('Continuous-time Fourier Transform');

运行效果如下：

如果想确认变换的正确性，可以在运行完上面这个脚本后，在命令行输入

plot(W/(2*pi*1000),(0.002./(1+(W./1000).^2))*1000);

xlabel('Frequency in KHz'); ylabel('Xa(jW)*1000');

运行效果如下：

这时会发现，根据上面推导的变换公式直接plot出的图形和变换后得到的图形是一样的，这样可以确定变换的正确性。

存在问题

目前存在的问题是，对于复函数的变换结果不正确。我想了很多天都找不出问题所在，只能暂时放弃，等以后有机会再研究。

离散时间傅里叶变换

下面是对上一个例子中的模拟输入信号做离散化，然后再进行离散傅里叶变换。

为了体现Nyquist定理，将使用两种不同的采样频率

使用Fs=5000sam/sec采样来获得x1(n)

使用Fs=1000sam/sec采样来获得x2(n)

% Analog Signal

Dt = 0.00005;

t = -0.005:Dt:0.005;

xa = exp(-1000*abs(t));

% Discrete-time Signal

Ts = 0.0002;

n = -25:1:25;

x = exp(-1000*abs(n*Ts));

% Discrete-time Fourier transform

K = 500;

k = 0:1:K;

w = pi*k/K;

X = x*exp(-j*n'*w); X = real(X);

w = [-fliplr(w), w(2:K+1)];

X = [fliplr(X), X(2:K+1)];

subplot(2,1,1);plot(t*1000,xa);

xlabel('t in msec.');

ylabel('x1(n)');

title('Discrete Signal');hold on;

stem(n*Ts*1000,real(x));gtext('Ts=0.2 msec');hold off;

subplot(2,1,2);plot(w/pi,X);

xlabel('Frequency in pi units');ylabel('X1(w)');

title('Discrete-time Fourier Transform');

Fs=5000sam/sec

xa(t)的频率为2KHz，因此它的Nyquist频率为4KHz，而它的采样频率为5KHz，所以是满足Nyquist采样定律的，此时不会发生混叠。

运行效果如下：

Fs=1000sam/sec

这里使用的采样频率为1KHz，不满足Nyquist条件，因此会发生混叠。观察一下就会发生，1KHz采样得到的序列的频域波形和前面的频域波形不同，这就是混叠导致的，而且过低的采样率采集的信号的变换的不可逆的。

运行效果如下：

收起

展开全文
【算法】FFT-快速傅里叶变换
2021-08-12 16:03:26
目录FFT作用数学知识代码：例题1：多项式乘法例题2：大数相乘碎碎念： FFT作用 FFT在算法竞赛中就有一个用途：加速多项式乘法多项式：形如 a0X0+a1X1+......+anXna_0X^0 +a_1X^1+......+a_nX^na0X0+a1X1+.........
目录

FFT作用
数学知识
代码：
例题1：多项式乘法
例题2：大数相乘

碎碎念：
FFT作用

FFT在算法竞赛中就有一个用途：加速多项式乘法

多项式：形如 $a_0X^0 +a_1X^1+......+a_nX^n$ 的代数表达式，可以记作 $f ( X ) =a_0X^0 +a_1X^1+......+a_nX^n$ ，其中， $a_0,a_1,a_2,......,a_n$ 是多项式的系数。

如果有两个多项式 $f (X), g (X)$ ，要求这两个多项式的乘积。那么最朴素的做法是每一位相乘，时间复杂度是 $O(n^2)$ ，很多情况下，需要优化，所以，就出现了FFT，时间复杂度是 $O (n l o g n)$

数学知识

多项式，复数，单位根，多项式的系数表达法，多项式的点值表达法…
有亿点点复杂，给一个链接，自己能理解多少就理解多少吧。
详解

代码：

其实FFT就是多项式的系数表示法与点值表示法之间的互换。
系数表达法容易用公式表示，但是相乘不容易；点值表示法相乘容易，不方便表示。
（点值表示法相乘：两个多项式P,Q分别取点 $x,y_1 )和( x, y_2 )$ ， $P * Q$ 就是点 $(x, y 1 * y 2)$ ,所以 $P * Q$ 的多项式的点值表示法就是 $(x, y 1 * y 2)$ ）
type=1：系数表示法 ——> 点值表示法
type=-1：点值表示法 ——> 系数表示法

例题1：多项式乘法

P3803 多项式乘法

题目：

代码：
```
#include <cstdio>
#include <cmath>
#include <iostream> 
using namespace std;
const int maxn=4*1e6+10;
const double pi=acos(-1.0);
struct Complex
{
    double x, y;
    Complex(double x=0,double y=0):x(x),y(y){}
    Complex operator+(const Complex &W) const
    {
        return {x + W.x, y + W.y};
    }
    Complex operator-(const Complex &W) const
    {
        return {x - W.x, y - W.y};
    }
    Complex operator*(const Complex &W) const
    {
        return {x * W.x - y * W.y, x * W.y + y * W.x};
    }
};
//a[i]表示当x=单位根的i次方时y的值 
Complex a[maxn],b[maxn];
int rev[maxn],n,m;
int l=1,bit=0;

void inif(int num)//l为位数,rev[i]代表i的位逆序置换操作后的结果 
{
	l=1,bit=0;
	while(l<=num) 
		l<<=1,bit++;
	for(int i=1;i<(1<<bit);i++)
		rev[i]=(rev[i>>1]>>1)|((1&i)<<bit-1);
}
void fft(Complex *f,int len,int type)
{
	for(int i=1;i<len;i++)
		if(rev[i]>i) swap(f[rev[i]],f[i]);
	for(int l=2;l<=len;l<<=1)//区间长度 
	{
		Complex wn=(Complex){cos(2*pi/l),type*sin(2*pi/l)};	//单位根 
		for(int i=0;i+l<=len;i+=l)
		{
			Complex w=(Complex){1,0};//幂 
			for(int k=i;k<i+(l>>1);k++,w=w*wn)
			{
				Complex t=w*f[k+(l>>1)],tmp=f[k];
				f[k]=tmp+t;//蝴蝶效应 
				f[k+(l>>1)]=tmp-t;
			}
		}
	}	
	if(type==-1)
	{
		for(int i=0;i<=n+m;i++) 
			f[i].x=(int)(f[i].x/l+0.5);
	}		 
}

int main ()
{
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&a[i].x);
	for(int i=0;i<=m;i++) scanf("%lf",&b[i].x);
	inif(n+m);
	fft(a,l,1); fft(b,l,1);//l是多项式项数
	for(int i=0;i<l;i++) a[i]=a[i]*b[i];
	fft(a,l,-1);
	for(int i=0;i<=n+m;i++) printf("%d ",(int)a[i].x);

	return 0;
}
```
例题2：大数相乘

P1919 A*B Problem升级版

题目：

解题思路：

对于每一个 n 位的十进制数，我们可以看做一个 n-1 次多项式 A，满足

$=a_0+a_1 \times 10+a_2\times10^2 +\cdots +a_{n-1}\times10^{n-1}$

那么对于两个大整数相乘,我们就可以卷起来辣！

小细节：
1. 首先要将读进去的大整数逆序存入多项式的系数中
2. 计算出相乘后的系数，别忘记处理进位问题
代码：
```
#include <cstdio>
#include <cmath>
#include <iostream> 
using namespace std;
const int maxn=4e6+10;
const double pi=acos(-1.0);
struct Complex
{
    double x, y;
    Complex(double x=0,double y=0):x(x),y(y){}
    Complex operator+(const Complex &W) const
    {
        return {x + W.x, y + W.y};
    }
    Complex operator-(const Complex &W) const
    {
        return {x - W.x, y - W.y};
    }
    Complex operator*(const Complex &W) const
    {
        return {x * W.x - y * W.y, x * W.y + y * W.x};
    }
};
Complex a[maxn],b[maxn];
int ans[maxn];
int rev[maxn],n,m;
int l=1,bit=0;

void inif(int num)
{
	l=1,bit=0;
	while(l<=num) 
		l<<=1,bit++;
	for(int i=0;i<(1<<bit);i++)
		rev[i]=(rev[i>>1]>>1)|((1&i)<<bit-1);
}
void fft(Complex *f,int len,int type)
{
	for(int i=0;i<len;i++)
		if(rev[i]>i) swap(f[rev[i]],f[i]);
	for(int l=2;l<=len;l<<=1) 
	{
		Complex wn=(Complex){cos(2*pi/l),type*sin(2*pi/l)};	
		for(int i=0;i+l<=len;i+=l)
		{
			Complex w=(Complex){1,0};
			for(int k=i;k<i+(l>>1);k++,w=w*wn)
			{
				Complex t=w*f[k+(l>>1)],tmp=f[k];
				f[k]=tmp+t;
				f[k+(l>>1)]=tmp-t;
			}
		}
	}	
	if(type==-1)
	{
		for(int i=0;i<=len;i++) 
			f[i].x=(int)(f[i].x/len+0.5);
	}		 
}


int main ()
{
	string s1,s2;
	cin>>s1>>s2;
	n=s1.size();
	n--;
	m=s2.size();
	m--;
	for(int i=0;i<=n;i++) //逆序存储 
		a[i].x=s1[n-i]-'0';
	for(int i=0;i<=m;i++) 
		b[i].x=s2[m-i]-'0';
	inif(n+m);
	fft(a,l,1); fft(b,l,1);
	for(int i=0;i<l;i++) a[i]=a[i]*b[i];
	fft(a,l,-1);
	for(int i=0;i<l;i++) //进位处理 
	{
		ans[i]+=a[i].x;
		if(ans[i]>=10)
		{
			ans[i+1]+=(ans[i]/10);
			ans[i]%=10;
		}
	}
	int pos=l;//从l开始找第一位不是0的数字 
	while(!ans[pos] && pos>0)	pos--;
	for(int i=pos;i>=0;i--)
		printf("%d",ans[i]);
	return 0;
}
```
碎碎念：

这个FFT，我只会用，但是对于背后神奇的数学推导，一知半解…
收起
展开全文
acm竞赛
【OpenCV 例程200篇】73. 二维连续傅里叶变换
2022-01-17 13:20:11

由于任何周期或非周期函数都可以表示为不同频率的正弦函数和余弦函数之和的形式，因此用傅里叶变换表示的函数特征完全可以通过傅里叶反变换来重建，而且不会丢失任何信息。这是频率域图像处理的数学基础。

收起

opencv python 图像处理计算机视觉
二维离散傅里叶变换手算实例（加OpenCV实现）
千次阅读 多人点赞 2020-11-06 19:53:57

二维离散傅里叶变换手算实例（加OpenCV实现）二维傅里叶原理及公式不做赘述，直接看图片理解吧如图想起来我还写了个计算二维离散傅里叶变换的代码要用opencv #include "opencv2/core.hpp" #include "opencv2/...

收起
二维离散傅里叶变换最详手算
千次阅读 2021-06-23 16:05:43

2 欧拉公式一定义二例题2.1 定义展开法2.2 矩阵累成法2.3 逆傅里叶变换法2.4 总结三参考文章前置知识 1. 1 复数定义(符号表示)： a + bi （这里a，b是实数，i是虚数单位，a叫做复数的实部，b叫做复数的...

收起
离散傅里叶变换（DFT）数学例子
万次阅读 多人点赞 2016-04-06 18:53:24

本篇博客主要是举个实例来展示离散傅里叶变换的计算过程（因为计算机主要是处理离散数值）

收起

离散傅里叶变换
使用matlab计算并绘制连续信号的傅里叶变换
千次阅读 2022-01-13 17:14:01

直接用定义式计算连续信号的傅里叶变换并绘图

收起

matlab
归纳傅里叶变换matlab代码_怎样用matlab进行傅里叶分析
2020-07-03 22:26:48

最新课件 %傅里叶变换 clc;clear all;close all; tic Fs=128%采样频率频谱图的最大频率 T=1/Fs%采样时间原始信号的时间间隔 L=256%原始信号的长度即原始离散信号的点数 t=(0:L-1*T%原始信号的时间取值范围 x=7*cos(2...

收起
【数字信号处理】傅里叶变换性质 ( 序列傅里叶变换共轭对称性质示例 )
千次阅读 2022-03-16 23:50:28

一、序列傅里叶变换共轭对称性质示例、 1、序列傅里叶变换共轭对称性质、 1、序列实部傅里叶变换、 2、序列虚部傅里叶变换、 3、共轭对称序列傅里叶变换、 4、共轭反对称序列傅里叶变换、 2、求 a^n u(n) 的傅里叶...

收起

傅里叶变换傅里叶变换性质共轭对称序列
快速傅里叶变换(fft + 例题 Gym100783C)
千次阅读 2020-05-17 17:49:12

快速傅里叶变换 (fast Fourier transform), 即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶...

收起
【数字信号处理】傅里叶变换性质 ( 傅里叶变换时移性质示例 )
千次阅读 2022-03-09 18:41:45

一、傅里叶变换线时移性质、二、傅里叶变换线时移性质示例

收起

傅里叶变换傅里叶变换时移性质
积分变换习题解答，包括傅里叶变换，拉普拉斯变换
2010-04-01 09:50:02

 dttdatfcossin02cos0所以 01||01||421021||22cossinttttfdt4、求下列函数的傅里叶变换。（1） 1||,01|||,|...

收起
第三章离散傅里叶变换及其快速算法习题答案参考.doc
2022-05-12 19:09:02

第三章离散傅里叶变换及其快速算法习题答案参考.doc

收起
傅里叶变换 二维离散傅里叶变换
万次阅读 热门讨论 2019-11-07 15:41:28

DFT:(Discrete Fourier Transform)离散傅里叶变换是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT的频域采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列...

收起

二维DFT
快速傅里叶变换（FFT）与快速数论变换（NTT）+例题
千次阅读 2018-08-07 16:09:49

因此，一种算法叫做“快速傅里叶变换”诞生了，它可以在 O ( n l o g n ) O ( n l o g n ) O(nlogn) 的时间内完成上述两部转化。快速傅里叶变换（FFT）单位根若 x n = 1 x n = 1 x^n=1 ，则 x x x ...

收起
二维傅里叶变换
2022-05-16 17:21:24

二维傅里叶变换 以下公式定义 m×n 矩阵 X 的离散傅里叶变换 Y。和是以下方程所定义的复单位根。 i 是虚数单位，p 和 j 是值范围从 0 到 m–1 的索引，q 和 k 是值范围从 0 到 n–1 的索引。在此公式中，X ...

收起

matlab
傅里叶变换
2022-03-11 13:45:27

傅里叶变换 傅里叶逆变换的傅里叶变换 线性性质 , 位移性质 1. 2. 微分性质如果在连续且只有有限多个可去间断点，且当时， ,则积分性质设,若则卷积的定义：已知两个函数，则称为的卷积。记作注意：卷积是一...

收起
【零基础】看懂理解傅里叶变换后的频谱图-附例题
万次阅读 多人点赞 2020-04-07 00:49:10

终于能看懂黑乎乎的傅里叶变换图了

收起

计算机视觉
离散时间傅里叶变换理论.pdf
2020-04-18 20:01:19

离散时间傅里叶变换,分7小节，全面完整推导DTFT公式，讨论其性质及其用于LTI系统分析的优势。课件例题丰富，内容完整，借鉴了诸多网站资源。参考书籍为欧本海姆的经典的信号与系统第二版。

收起
2005_09傅里叶变换及其应用（第3版）_11482158.pdf
2021-05-08 06:55:46

2005_09傅里叶变换及其应用 (第3版)_11482158[General Information]书名=2005.09傅里叶变换及其应用 (第3版)作者=(美)罗纳德·N·布雷斯韦尔著页数=480SS号出版日期=2005.09前言目录译者序作者简介前言第1章　绪论第...

收起
常用的傅里叶变换 定理各种变换的规律(推荐).pdf
2019-09-15 19:48:23

常用的傅里叶变换 定理各种变换的规律(推荐)

收起
对二维图像进行离散傅里叶变换
千次阅读 2022-07-07 19:46:27

本文基于二维离散傅里叶变换，实现对灰度图像的离散傅里叶变换，包括有Matlab自带的函数，以及自己编写的函数。

收起

python
《数字信号处理》——离散傅里叶变换与z变换专题
千次阅读 2021-04-23 21:05:21

离散非周期信号的傅里叶变换（DTFT）离散时间傅里叶变换(DTFT，频率连续) 离散傅里叶变换(DFT) 快速傅里叶变换(FFT，DFT的优化) 基本概念连续信号与离散信号从横轴时间上分析模拟信号与数字信号从纵轴幅值上分析 ...

收起
二维离散傅里叶变换
千次阅读 2021-11-19 10:39:43

二维离散傅里叶变换的直观理解一维离散傅里叶变换很简单，就是将一个向量（一维信号）分解为多个频率的向量（一维信号之和）。那么二维离散傅里叶变换我们类比一下就知道了，他只是将一幅图片（二维信号）分解为多...

收起

人工智能
第三章 傅里叶变换和离散傅里叶级数（DFS）
2021-10-04 11:38:06

关于傅里叶变换和离散傅里叶级数求和范围的理解理解： important！！！ 傅里叶的精髓是基函数的正交展开，用n项独立不相关的基频函数可以线性组合出目标函数就够啦！基于这一思想，来理解为何有时候积分范围是(−...

收起