-
C语言快速傅里叶变化FFT
2020-04-19 14:59:56快速傅里叶变化:主要有几个难点 1.如果使用DIT首先要倒位序,网上的方法值得参考 思路是和二进制进位的思路相同只是是反着进位 2.要自己封装复数类型 3.难点是公式里的k不好求,因为k不是连续取得的,我用了一个数j...快速傅里叶变化:主要有几个难点
1.如果使用DIT首先要倒位序,网上的方法值得参考
思路是和二进制进位的思路相同只是是反着进位
2.要自己封装复数类型
3.难点是公式里的k不好求,因为k不是连续取得的,我用了一个数j来存储,这样可以把k计算出来
4.wn^r中r的取值,直接移位即可,不需要转化为二进制
5.输入不是2的倍数时自动补零
其余我都写在注释里了,感觉还是有很多地方要改进#include <stdio.h> #include <math.h> #include <stdlib.h> #include "复数类型.c" //可以用M_E表示自然对数,M_PI表示圆周率 /*首先要导位序,也可以只要一个存储空间A[i],只需要改变A[i]中i 重点1.只改变当i小于j的时候的值,而i大于等于j时因为已经换过,不能再换回来了 重点2.思路是与二进制算法思路,小于N/2时即,最高位为0,此时的结果需要加上N/2 当大于N/2时,最高位为1,此时要向低位进位,意味着要先-N/2代表进位 再判断是否大于N/2,大于再向次低位进位,否则直接加N/4*/ void daoxu(double *x,int n) { int i=0,j=0,k=0;//i表示原始值,而j表示倒位序后的值 double temp; for(;i < n;i++) { if(i<j) { temp = x[j]; x[j] = x[i]; x[i] = temp; } k = n/2; while(j >= k && k != 0) { j = j - k; k = k/2; } j = j + k; } /*for(i=0;i<n;i++) { printf("x[%d]=%f\n",i,x[i]); //printf("result: new_k[%d]=%d\n",i,new_k[i]); }*/ } int process(int* ori_len)//对输入的数据进行处理,还要对数据进行补零,问题在于如何分配内存,使用链表? { int temp = 1,i; for(i=0;;i++) { if(temp >= *ori_len) break; temp =temp * 2; } *ori_len = temp; return i;//返回幂次方 } complex w_n(int k,int m,int L,int length) //如果这个值是复数意味着要包装复数类型,并且和这个数做四则运算的其他数也必须是复数 //需要新建复数这个数据结构,但是其实最后计算的是幅值 { int r=k<<(L-m);//不应该取这个k,而应该取倒位序的k r=r&(length-1);//去除高位 //printf("得到的r为=%d\n",r); complex result; result.real=cos(2*M_PI/length*r);//用欧拉公式展开 //printf("result.real=cos(2*pi/%d*%d)=%f\n",length,r,result.real); result.imag=-sin(2*M_PI/length*r); //printf("result.imag=-sin(2*pi/%d*%d)=%f\n",length,r,result.imag); return result; } void butterfly_operation(double *x,int powe,int length) { int k_num = length/2,k_cal,k,l; complex r; complex *back_up=(complex*)calloc(length,sizeof(complex));//需要备份值来保留处理前的数据,定义为虚数类型 complex *y=(complex*)calloc(length,sizeof(complex)); for(int m=1;m<=powe;m++)//pow表示蝶形的次数 { //printf("第%d次蝶形运算\n",m); //printf("要做%d次蝶形运算\n",k_num); k = 0;//运算之前将k清零,这里的k是实际是倒位序的 l = 0; for(int i=0;i<length;i++)//运算完一次蝶形更新备份值 { if(m==1) back_up[i] = cpx_make(x[i],0); else back_up[i] = y[i]; //printf("back_up[i]=%f+%fi\n",back_up[i].real,back_up[i].imag); } k_cal=pow(2,m-1);//k_cal实际代表的是2^(m-1)的值 //printf("计算2^(m-1)的值=%d\n",k_cal); for(int j=0;j<k_num;j++)//在做第m次蝶形的时候,要做k_num次蝶形运算 { //printf("k=%d\n",k); r = w_n(k,m,powe,length); //printf("r=%f+%fi\n",r.real,r.imag); y[k] =cpx_add(back_up[k],cpx_mul(back_up[k+k_cal],r));//x[k]取的是模值 //printf("y[%d]=%f+%fi\n",k,y[k].real,y[k].imag); y[k+k_cal] = cpx_sub(back_up[k],cpx_mul(back_up[k+k_cal],r)); //printf("cpx_mul(back_up[k+k_cal],r)=%f+%fi\n",cpx_mul(back_up[k+k_cal],r).real,cpx_mul(back_up[k+k_cal],r).imag); //printf("y[%d]=%f+%fi\n",k+k_cal,y[k+k_cal].real,y[k_cal+k].imag); if(k_cal>=2&&k<l+k_cal-1) { k+=1; } else { k=k+k_cal+1; l=k; } } //printf("-------------------------------------------------------------------\n"); } for(int i=0;i<length;i++) { x[i]=sqrt(y[i].real*y[i].real+y[i].imag*y[i].imag); //printf("y[i].real=%f,y[i].imag=%f\n",y[i].real,y[i].imag); } free(back_up); free(y); } int main() { int length,power,old_length;//以后要更正 double t; printf("请输入采样的长度length="); scanf("%d",&length);//输入采样的长度*/ double *x,*y; y=(double*)calloc(length,sizeof(double));//可以动态分配内存 for(int i=0;i<length;i++) { //y[i] = i+1;//在输入之前就会分配好内存空间因此不能临时输入 t=i*0.001; y[i] = 1.2+2.7*cos(2*M_PI*33*t)+5*cos(2*M_PI*200*t+M_PI/2); //y[i]=0.6*sin(2*M_PI*500*i)+0.6*sin(2*M_PI*50*i); //y[i] = sin(M_PI*i/8); //printf("%f\n",y[i]); } old_length = length; power = process(&length);//改变了length x=(double*)calloc(length,sizeof(double)); for(int i=0;i<old_length;i++) { x[i] = y[i]; printf("x[i]=%f\n",x[i]); } for(int i=old_length;i<length;i++) { x[i]=0; printf("x[i]=%f\n",x[i]); } free(y); daoxu(x,length); butterfly_operation(x,power,length); for(int i=0;i<length;i++) printf("%d\t%f\n",i,x[i]);//在输入之前就会分配好内存空间因此不能临时输入 free(x); return 0; }
我自己封装的复数头文件
#include <stdio.h> #define Real(c) (c).real #define Imag(c) (c).imag typedef struct { double real; double imag; }complex; complex cpx_make(double real, double imag) { complex ret; ret.real = real; ret.imag = imag; return ret; } complex cpx_add(complex a,complex b) { return cpx_make(Real(a) + Real(b), Imag(a) + Imag(b)); } complex cpx_sub(complex a,complex b) { return cpx_make(Real(a) - Real(b), Imag(a) - Imag(b)); } complex cpx_mul(complex a,complex b) { return cpx_make(Real(a) * Real(b)-Imag(a) * Imag(b),Real(a) * Imag(b) + Imag(a) * Real(b)); } complex cpx_div(complex a,complex b) { Real(b) = 1/Real(b); Imag(b) = 1/Imag(b); return cpx_mul(a,b); }
-
《数据科学家养成手册》傅里叶变换与反傅里叶变换笔记
2018-04-28 18:25:52最早的电话使用的模拟信号原理图(1)...局限性就是当距离比较远 的时候就很难还原(电流的衰减大)傅里叶变换和反傅里叶变换就解决了这个问题--使得数字信号通信代替模拟信号通信称为可能(1)任意一个是与信号f(...最早的电话使用的模拟信号原理图
(1)声音通过金属振动膜感应声波来影响磁场和电流,并将这种带有金属振动膜振动的“信息”的电流传递给另一端
(2)另一端则进行反向工作,把不断变化的电流转化为电线圈中磁场的变化,使得金属振动膜发出同样的振动。
局限性就是当距离比较远 的时候就很难还原(电流的衰减大)
傅里叶变换和反傅里叶变换就解决了这个问题--使得数字信号通信代替模拟信号通信称为可能
(1)任意一个是与信号f(t)经过傅里叶变换,都可以转化为多个余弦波叠加的形式
一旦这个环节被我们掌控的话,语音信号就可以通过调制,由高频载波用极高的速率发送电信号和光信号。
而在接收端,用滤波器过滤余弦波信号后,使用傅里叶反变换进行是与信号的还原
对傅里叶的理解资源:
深入浅出的讲解傅里叶变换(真正的通俗易懂) - CSDN博客
https://blog.csdn.net/l494926429/article/details/51818012实验的来源:
实验六傅里叶变换及其反变换_百度文库
https://wenku.baidu.com/view/1ba79654178884868762caaedd3383c4bb4cb4bd.html代码:
syms t v w x phase im re ; %定义符号变量 a = input('请输入a=') f = exp(-a*abs(t)); %f(t) = exp(-2*t)*u(t) Fw = fourier(f); %求傅里叶变换 subplot(311); ezplot(f); %绘制f(t)的时域波形 axis([-1 2.5 0 1.1]); subplot(312); ezplot(abs(Fw)); %绘制幅度谱 im = imag(Fw); %计算F(w)的虚部 re = real(Fw); %计算F(w)的实部 phase = atan(im/re); %计算相位谱 subplot(313); ezplot(phase); %绘制相位谱
syms t v w x phase im re ; %定义符号变量 F(w) = 1/(1+w^2); %待求解变换的公式 f = fourier(Fw,t); subplot(311); ezplot(f); %绘制f(t)的时域波形 axis([-1 2.5 0 1.1]); subplot(312); ezplot(abs(Fw)); %绘制幅度谱 im = imag(Fw); %计算F(w)的虚部 re = real(Fw); %计算F(w)的实部 phase = atan(im/re); %计算相位谱 subplot(313); ezplot(phase); %绘制相位谱
-
快速傅里叶变换
2016-11-30 14:27:13快速傅里叶变化(FFT)是离散傅里叶变化(DFT)的快速算法。一个长度为N的有限长序列的DFT为 式中的 为 。离散傅里叶反变换为 FFT算法的基本思想是利用 的对称性和周期性,即 =- 和 (r为任意整数),通过将长...快速傅里叶变化(FFT)是离散傅里叶变化(DFT)的快速算法。一个长度为N的有限长序列的DFT为
式中的 为 。离散傅里叶反变换为
FFT算法的基本思想是利用 的对称性和周期性,即 =- 和 (r为任意整数),通过将长序列分解为短序列,再由短序列的DFT组合来实现整个长序列的DFT,从而减少DFT的复数乘法和复数加法次数,提高DFT的运算速度和效率。将N点的长序列分解为短序列可以有多种方式,常用的有基2时间抽取FFT算法,基2频率抽取FFT算法,基4抽取FFT算法和分裂基FFT算法。基2抽取算法简单有效,基4抽取算法和分裂基抽取算法与基2抽取算法原理类似,但基4算法和分裂基算法效率更高,运算量更少。本文主要介绍基2抽取FFT算法。
1、基2时间抽取FFT算法:
基2时间抽取FFT算法是将N点序列x[n]按序号奇偶分解得到两个 点子序列 和 ,通过计算这两个子序列的DFT 和 实现 的DFT 。算法要点为N点序列x[n]对时域下标按奇数和偶数分解,即
(1)
对频域X[k]按前一半和后一半分解,即 。频域的前一半可由下式求得
由DFT的周期性和 的对称性有
-=
故频域的后一半可由下式求得
将前一半和后一半组合起来,即得原来的长序列x[n]的DFT 。
可对 和 继续按式(1)分解下去,直到每个子序列的长度为1为止。任何一个N= 点的DFT,可由通过M次分解,得到长度为1的子序列。对于N不是2的整数次幂的情况,可在序列后面添0补齐到2的整数次幂。每级分解运算可用一个蝶形运算流程图来表示,如图
下图是对一个长度为8的序列进行基2时间抽取FFT运算的流程图
直接计算N点序列的DFT需要N2次复数乘法,需要N(N-1)次复数加法。而用基2时间抽取FFT算法,则乘法次数为 ,加法次数 ,运算次数大大减少。
2、基2频率抽取FFT算法
基2频率抽取FFT算法则是对频域下标按奇数和偶数分解,对时域下标按前一半和后一半分解。 的DFT 可表示为
将 分解成 和 两个长度为 的子序列,如下:
(2)
则 的DFT 可按频率下标分奇数和偶数表示为
可继续将 和 按式(2)分解下去,直到每个子序列长度为1。任何一个N= 点的DFT,可由通过M次分解,得到长度为1的子序列。对于N不是2的整数次幂的情况,可在序列后面添0补齐到2的整数次幂。每级分解运算也可用一个蝶形运算流程图来表示,如图
下图是对一个长度为8的序列进行基2频率抽取FFT运算的流程图:
3、编程实现
Matlab有一个内置的函数fft()可以实现快速傅里叶变换,由于它是built-in函数,所以无法查看源码来判断它是应用的那种分解方式来实现的快速傅里叶变换。下面将用Matlab语言来编写实现类似fft()的快速傅里叶变换函数myfft()。本函数采用的FFT算法为基2时间抽取FFT算法,函数定义源码如下:
%----------------------------------
%myfft(A,M)有两个形参,A为要进行快速傅里叶变换的序列,M为序列长度对2的对数,
%即如果序列A[n]的长度为8,则M=3
function [A] = myfft(A,M)
N=2^M; %序列的长度为N
LH=N/2; %序列长度的一半
N1=N-2;
J=LH
for I=1:1:N1 %将输入序列A[n]按运算流程图第一列的顺序排好
if I<J
temp=A(I+1);
A(I+1)=A(J+1);
A(J+1)=temp;
end
K=LH;
while J>=K
J=J-K;
K=K/2;
end
J=J+K;
end
for L=1:1:M %L表示序列分解层次,此次分解有M-L个子序列
B=2^(L-1);
for S=0:B-1 %执行第S个蝶形运算,每个子序列共有B个蝶形运算
p=S*2^(M-L); %旋转因子的上标,即此时旋转因子为
for k=S:2^L:N-1 %执行各个子序列中的第S个蝶形运算,有M-L个子序列
temp=A(k+1)+A(k+B+1)*exp(-i*2*pi*p/N);
A(k+B+1)=A(k+1)-A(k+B+1)*exp(-i*2*pi*p/N);
A(k+1)=temp; %蝶形运算
end
end
end
end
python语法简洁清晰,具有强大和丰富的三方库,能把其他语言(如C/C++,Matlab,java等)制作的各种模块轻松联结在一起,开发效率高,很适合快速开发实验原型以验证算法正确性和有效性。下面是基2时间抽取FFT算法的一个python版本实现。
#-------------------------------------
import numpy as np
def FFT(x):
x =np.asarray(x, dtype=float)
N =x.shape[0]
ifN%2>0: #输入序列的长度必须为2的幂
raiseValueError("size of x must be a power of 2")
elif N<= 2: # 子序列长度小于等于2时,直接计算
x =np.asarray(x, dtype=float)
N = x.shape[0]
n = np.arange(N)
k = n.reshape((N, 1))
M = np.exp(-2j * np.pi * k * n / N)
eturn np.dot(M, x)
else:
X_even= FFT(x[::2]) #偶序列继续分解
X_odd= FFT(x[1::2]) #奇序列继续分解
factor= np.exp(-2j * np.pi * np.arange(N/2) / N) #旋转因子
returnnp.concatenate([X_even + factor* X_odd, #返回序列频域前半部分
X_even - factor*X_odd]) #返回序列频域后半部分
4、函数测试和实验
快速傅里叶变换可以得到信号的频谱。通过对一个给定的信号进行快速傅里叶变换,观察它的频谱,并和Matlab内置函数fft()变换得到的结果进行对比,来对myfft()的功能进行测试。测试可用如下脚本来实现:
%----------------------------
%测试脚本,分别用myfft()和fft()对一个信号进行快速傅里叶变换,该信号由一个
%50Hz和120Hz的正弦信号叠加,频谱会在这50Hz和120Hz出现峰值。
Fs = 1000; % 采样频率
T = 1/Fs; % 采样周期
L = 5000; % 信号长度
t = (0:L-1)*T; % 时间序号数组
% x(t)=0.7sin(2 *50t)+sin(2 *150t),为50Hz和150Hz正弦信号的叠加
x = 0.7*sin(2*pi*50*t) + sin(2*pi*150*t);
y = x + 2*randn(size(t)); % 给信号加上噪声
%画出x(t)和y(t)的图像
subplot(2,1,1)
plot(Fs*t(1:50),x(1:50))
title('Sinusoids Signal')
xlabel('time (milliseconds)')
subplot(2,1,2)
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
%对加了噪声的y(t)进行快速傅里叶变换
NFFT = 2^nextpow2(L); %延长信号长度为最小的2的幂
y=[y zeros(1,2^nextpow2(L)-L)]; %延长信号长度,后面用0补齐
Y = myfft(y,log2(NFFT))/L; % 用自己实现的快速傅里叶变化myfft处理信号得到Y(f)
Y2 = fft(y,NFFT)/L; %用Matlab官方实现的fft处理信号得到Y2(f)
f = Fs/2*linspace(0,1,NFFT/2+1);%频域横坐标
%对信号y(t)进行快速傅里叶变换,并画出的频谱
figure
subplot(2,1,1)
plot(f,2*abs(Y(1:NFFT/2+1))) %
title('myfft()处理的y(t)频谱')% y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
subplot(2,1,2)
plot(f,2*abs(Y2(1:NFFT/2+1))) %
title('fft()处理的y(t)频谱')
xlabel('Frequency (Hz)')
ylabel('|Y2(f)|')
%对原信号x(t)进行快速傅里叶变换,并画出频谱
x=[x zeros(1,2^nextpow2(L)-L)]; %延长信号长度,后面用0补齐
X = myfft(x,log2(NFFT))/L; % 用自己实现的快速傅里叶变化myfft处理信号得到X(f)
X2 = fft(x,NFFT)/L; %用Matlab官方实现的fft处理信号得到X2(f)
figure
subplot(2,1,1)
plot(f,2*abs(X(1:NFFT/2+1))) %
title('myfft()处理的x(t)频谱')%
xlabel('Frequency (Hz)')
ylabel('|X(f)|')
subplot(2,1,2)
plot(f,2*abs(X2(1:NFFT/2+1))) %
title('fft()处理的x(t)频谱')
xlabel('Frequency (Hz)')
ylabel('|X2(f)|')
运行脚本,得到x(t)和y(t)的频谱,如图。可见,myfft()和Matlab内置函数fft()的处理结果一致,没有区别。频谱图在50Hz和150Hz出现峰值,幅值分别为0.63和1.0左右。
Python版的测试脚本如下
a=np.array([1,0,1,0,1,0,1,0]) #定义序列a[n]=[1,0,1,0,1,0,1,0]
b=[1,1,1,1,0,0,0,0] #定义序列b[n]=[1,1,1,1,0,0,0,0]
A=FFT(b) #对a[n]进行快速傅里叶变换
B=FFT(a) #对b[n]进行快速傅里叶变换
printA,”\n”,B #输出两个序列的傅里叶变换的结果A[M]和B[m]
python版的傅里叶变换结果和Matlab自带的fft()变换结果对比如下
Python版的运行结果
Matlab版的运行结果
注意python版本的运行结果中,有”e-16”的结果,这是因为计算机浮点数的运算误差导致的,实际上可认为是0。可见,python实现的快速傅里叶变换函数的结果和Matlab的fft()变换结果一致,python实现的快速傅里叶变换结果正确无误。
5、快速傅里叶变换的应用
快速傅里叶变换可以将有限长或者周期离散时间信号从时域表示转换到频域表示。快速傅里叶变换除了可以得到一个信号的频谱外,还可以基于得到的频谱对信号进行一些频域上的处理,并通过快速傅里叶反变换得到经过处理后的信号。比如计算离散时间信号的线性卷积,滤波,信号压缩,降噪,谱估计等。FFT除了用于离散时间信号处理,也可以用于逼近连续时间信号的频谱。
下面脚本展示了利用快速傅里叶变换和反变换进行快速卷积计算
subplot(3,1,1);
n=0:16;
x=0.8.^n; %x[n]=0.8^n的离散信号
stem(n,x);
xlabel('n');ylabel('x[n]');
subplot(3,1,2);
n=0:15;
y=[ones(1,10) zeros(1,6)];
stem(n,y);
xlabel('n');ylabel('y[n]')
subplot(3,1,3);
L=26;n=0:L-1;
X=fft(x,L);
Y=fft(y,L);
Z=X.*Y; %在频域两信号的频谱相乘
z=ifft(Z,L); %快速傅里叶反变换得到两信号的卷积
stem(n,z);
xlabel('n');
ylabel('z[n]')
脚本运行结果如图
在工程领域,FFT可以用于语音处理,工程传感器信号处理,特征提取等。
二维FFT和IFFT可应用到图像处理领域,比如图像降噪,模糊化,清晰化,轮廓化,压缩,还原等。快速傅里叶变换及其反变换是图像处理领域的重要处理手段,是进行其他各项处理的基础,应用场合很广。下面脚本展示了图像的二维快速傅里叶变换及其反变换,及基于快速傅里叶变换得到的频谱图,相位图等一系列频域信息。
C=imread('2.jpg');%载入图片
A=rgb2gray(C);
B=fftshift(fft2(A));% 进行傅立叶变换
E=fft2(A);
D=ifft2(E);
subplot(331)
imshow(C);
title('原始图像');
subplot(332)
imshow(A);
title('原始灰度图像');
subplot(339)
imshow(abs(B),[]);
title('原始频谱图');
subplot(334)
imshow(log(abs(B)),[]);
title('取对数后的频谱图');
subplot(335)
imshow(angle(B),[]);
title('相位图');
subplot(336)
imshow(real(B),[]);
title('实部图');
subplot(337)
imshow(imag(B),[]);
title('虚部图');
subplot(338)
imshow(abs(E),[]);
title('fft2频谱图');
subplot(333)
imshow(D,[]);
title('经快速傅里叶逆变换得到的灰度图');
下面为脚本运行的结果,可见经快速傅里叶反变换后,可以从图像的频域表示得到时域表示,即重新得到一张图像。由于此脚本没有对图像做其他处理,所以快速傅里叶反变换得到的图形与原来的灰度图像一致。
6、总结
快速傅里叶变换大大简化了运算量,并且算法具有迭代和分治的特点,特别适合计算机实现,使得用计算机对离散时间信号进行快速处理成为现实。配合采样定理和信号采样及重构等方法,使得快速傅里叶变换可以进一步推广到连续时间信号处理,大大提高了信号处理的效率和速度。在语音处理领域,傅里叶变换可以用于语音识别,音频分离,降噪等应用。在图像应用领域,快速傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数,傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。而灰度分布函数即为灰度图的表达。快速傅里叶变换使得对图像可以用频域处理的方法来处理,对基于计算机的图像处理应用的可行性和实时性具有重大意义,是数字信号处理的基石。
-
图像傅里叶变换的幅度谱和相位谱的以及反变换
2017-05-28 16:27:42目的:读取图像 A(lena.tiff)和B(rice.tif),显示这两幅图像,对图像作傅立叶变换,显示图像的傅里叶幅度谱和相位谱。做傅立叶逆变换,显示重建图像。图像的频率是表征图像中灰度变化剧烈程度的指标,是灰度在平面...目的:读取图像 A(lena.tiff)和B(rice.tif),显示这两幅图像,对图像作傅立叶变换,显示图像的傅里叶幅度谱和相位谱。做傅立叶逆变换,显示重建图像。
图像的频率是表征图像中灰度变化剧烈程度的指标,是灰度在平面空间上的梯度
对图像而言,图像的边缘部分是突变部分,变化较快,因此反应在频域上是高频分量;图像的噪声大部分情况下是高频部分;图像平缓变化部分则为低频分量。也就是说,傅立叶变换提供另外一个角度来观察图像,可以将图像从灰度分布转化到频率分布上来观察图像的特征。
imshow()函数:
imshwo()函数用于接收一个像素矩阵,显示该图像,其显示的参数有两种类型
unit8;像素在矩阵处理范围为0-255
double:若值大于1,转化为1,若小于1,转化为0图像进行二维傅立叶变换得到频谱图,就是图像梯度的分布图,当然频谱图上的各点与图像上各点并不存在一一对应的关系,即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点,实际是图像上某一点与邻域点差异的强弱,即梯度的大小,也即该点的频率的大小(可以这么理解,图像中的低频部分指低梯度的点,高频部分相反)。
代码如下:%%图像的傅里叶变换%% imA=imread('rice.tif','tif'); %读取图像 imB=imread('lena.tiff','tif'); subplot(2,3,1); imshow(imA); title('原图像A'); subplot(2,3,2); imshow(imB); title('原图像B'); FA=fft2(imA);%对图像进行傅里叶变换 FB=fft2(imB); fA=fftshift(FA); %对图像频谱进行移动,是0频率点在中心 fB=fftshift(FB); sA=log(abs(fA));%获得傅里叶变换的幅度谱 sB=log(abs(fB)); phA=log(angle(fA)*180/pi);%获得傅里叶变换的相位谱 phB=log(angle(fB)*180/pi); subplot(2,3,3); imshow(sA,[]); %显示图像的度谱,参数与[]是为了将sA的值线形拉伸 title('图像A的傅里叶变换幅度谱'); subplot(2,3,4); imshow(phA,[]); %显示图像傅里叶变换的相位谱 title('图像A傅里叶变换的相位谱'); subplot(2,3,5); imshow(sB,[]) title('图像B的傅里叶变换幅度谱'); subplot(2,3,6); imshow(phB,[]); title('图像B傅里叶变换的相位谱'); A=ifft2(FA);%傅里叶反变换 B=ifft2(FB); figure subplot(1,2,1); imshow(A,[]); title('傅里叶反变换得到的A图像'); subplot(1,2,2); imshow(B,[]); title('傅里叶反变换的到的B图像');
结果如下:
-
案例解释图像傅里叶变换的幅度谱和相位谱的以及反变换
2017-11-09 13:21:01目的:读取图像 A(lena.tiff)和B(rice.tif),显示这两幅图像,对图像作傅立叶变换,显示图像的傅里叶幅度谱和相位谱。做傅立叶逆变换,显示重建图像。 图像的频率是表征图像中灰度变化剧烈程度的指标,是灰度在... -
numpy 傅里叶变换与反变换高低通滤波与带通滤波
2017-04-17 14:53:06#coding=utf-8 import cv2 import numpy as np import matplotlib.pyplot as plt img=cv2.imread('test1-angle.jpg',cv2.IMREAD_GRAYSCALE) # f = np.fft.fft2(img) ...# #取绝对值:将复数变化成 -
有关傅里叶变换知识
2021-03-16 12:12:25傅里叶变换是一种函数在空间域和频率域的变化,从空间域到频率域的变换时傅里叶变换,而从频率域到空间域的变换是傅里叶的反变化。 频域:是指在对函数或信号进行分析时,分析其和频率有关 部分,而不是和时间有关... -
图像的傅里叶变换
2020-07-21 17:41:24傅里叶变换及其反变换是线性处理 低频——图像中变化缓慢的灰度分量 高频——图像中变化急剧的灰度分量,例如边缘和噪声 频率域滤波过程:修改图像的傅里叶变换,再计算IDFT返回到图像域 低通滤波器——衰减高频... -
基于傅里叶变换图像配准
2019-01-07 11:05:22如果图像是图像经平移后的图像,即,则对应的傅里叶变化F1和F2的关系为: (1) 且对应频域中两个图像的互能量谱为: (2) 式中为的复共轭。平移理论表明,互能量谱的相位等于图像间的相位差。通过对互能量谱进行反... -
常用傅里叶变换公式表_第四章 频率域滤波-(二)基本概念之连续变量函数的傅里叶变换...
2020-12-11 08:32:22我们啰嗦了这么久,从复数讲到欧拉公式,从级数讲到冲激函数,今天终于可以一睹傅里叶变化的风采了。连续函数傅里叶变换描述总是空洞的,直接上公式。 由 表示的连续变量 的连续函数 的傅里叶变换公式。其中 也是一... -
傅里叶变换原理---OpenCV-Python开发指南(31)
2021-02-23 17:51:20目录前言傅里叶变化实现傅里叶变化实现逆傅里叶变化高通滤波与低通滤波 前言 要理解傅里叶变换,我们首先需要了解图像处理。在图像处理的过程中,一般分为空间域处理和频率域处理。 空间域处理是直接对图像内的像素... -
求周期方波信号的傅里叶级数_傅里叶变换
2020-12-31 23:05:08变换为了能说傅立叶变换,我们从一个简单的变换开始说起在平面直角坐标系中有一个向量如图:那么我们可以用(4,3)来表示a向量,用(1,3)来表示b向量,这样我们把图变化成了一个数,反过来,我们可以用一个数表示一个图,... -
Halcon:傅里叶变换(4)高通低通滤波器
2021-01-11 13:46:44反傅里叶变换后:边缘,菱角就不清晰了,类似模糊了一样 高通滤波器:代表高频通过的滤波器,就是中间的低频不通过,四周的高频通过。 正常表现形式就是如下: 中间为黑色通过区域,背景为白色不通过区域 反傅里叶... -
matlab 通过矩阵变换使图像旋转平移_光电图像处理 | 傅里叶变换(一)
2020-11-28 10:53:21傅里叶变换Fourier transform 1 傅里叶变化基本知识1.1 一维连续Fourier变换对函数f(x)进行傅里叶变换得到F(u)逆变换:从F(u)到f(x)进行反傅里叶变换一维连续函数f(x)的傅立叶变换F(u) 一般是虚数,可用复数形式表示... -
有关数字图像处理 图像的傅里叶变换
2017-11-16 08:35:061、将一幅图分别进行X轴与Y轴上的平移,所得的傅里叶频谱与原图像频谱的傅里叶频谱有什么变化,请说明理由。 2、将一幅图进行离散傅里叶变换,得到其傅里叶频谱图,在对原图像进行一定角度的旋转,得到的频谱图与原... -
opencv傅里叶变换图片显示问题
2021-01-19 22:00:56# 第三步 使用cv2.dft进行傅里叶变化 dft_image = cv2.dft(float_image, flags=cv2.DFT_COMPLEX_OUTPUT) # 第四步 使用np.fft.fftshift 将低频部分转换到图像的中心 shift_image = np.fft.fftshift(dft_... -
-
数字图像处理,图像的频域变换(二)——二维离散傅立叶变换的性质。都是数学都是各种傅里叶各种性质
2019-07-24 21:35:33https://zh.wikipedia.org/wiki/%E6%AC%A7%E6%8B%89%E5%85%AC%E5%BC%8F ...可以先对y轴进行傅里叶变化,然后对x轴处理。反过来也一样 获得的图... -
下一代无线局域网(802.11n) 第2章 OFDM
2020-10-13 16:13:31一个基带OFDM波形是由一组系数XkX_kXk的反傅里叶变化组成的,即 r(t)=1N∑kXkej2πkΔFt 0<t<T r(t) = \frac{1}{N} \sum_{k} X_k e^\mathrm{j2\pi k \Delta_F t} ~~ 0<t<T r(t)=N1k∑... -
计算机图形学笔记(五):反走样
2020-08-20 12:50:46采样带来的问题 锯齿、摩尔纹、Wagon轮效应 ...只留下高频信息,再反傅里叶变换,得到的图像主要是原图像中物体的边界 低通滤波 只留下低频信息,再反傅里叶变换,得到的图像是原图像的模糊 更多滤 -
算法在ros中应用_原力科普之 人工智能算法在超导测量中的应用
2020-12-17 05:00:44毫无疑问人工智能、深度学习是现在非常火热的学科,会是报考的热点。...该算法在图像识别领域大量应用,其核心思路是利用傅里叶变化将时域中的卷积与反卷积变化为频域中的乘法与除法,对时域中的看似... -
7/13数字图像处理笔记
2020-07-13 20:05:507月12日频率域滤波基本概念傅里叶级数傅里叶变换冲激及其取样特征卷积 频率域滤波 滤波器:抑制或最小化某些频率的波或振荡的装置或材料 频率:自变量单位变化期间,一个周期函数重复相同值序列的次数 ...傅里叶反变 -
mysql 把 case的结果作为变量引用_Maxwell中FFT结果作为输出变量方法
2020-12-19 00:39:06使用软件的FFT功能可以快速将瞬态分析的结果曲线进行频谱输出。在电机设计过程中需要对反...一、原理分析傅里叶变化的余弦分量与正弦分量分别为: 则 次分量幅值为: 由于软件中提供了avg函数,avg(f(t))函数定义... -
实验五 Wav信号的波形分析.doc
2020-04-24 22:18:53PAGE PAGE 1 实验五 Wav信号的波形分析 一 实验目的 借助本实验帮助同学们巩固傅里叶变换及其反变换的知识学习从时域和频域两个角度来观察信号并尝试利用短时傅里叶变换分析非平稳信号的频谱变化 二 实验原理 借助... -
补充知识——数字图像处理(3)
2021-01-30 21:13:49傅里叶变换是一种函数在空间域和频率域的变换,从空间域到频率域的变换是傅里叶变换,而从频率域到空间域是傅里叶的反变换。 时域与频域: 频域(frequency domain):是指在对函数或信号进行分析时,分析其和频率... -
图神经网络(七)A Generalization of Convolutional Neural Networks to Graph-Structured Data
2021-03-13 19:59:16在谱域卷积中,我们的基本思路是:先将空域输入信号和空域卷积核转换到谱域,然后在谱域中相乘,再通过傅里叶反变换转换回空域。 图结构中信号的空域和谱域的转换是由图傅里叶变换实现的,其中拉普拉斯特征向量作为... -
数字图像处理复习(part2)
2019-12-25 12:07:20一维离散傅里叶变换及其反变换 模拟图像和数字图像 模拟图像:又称连续图像,是指在二维坐标系中连续变化的图像,图像的像点是无限稠密的,同时具有灰度值(明暗变化)。 Chapter1 绪论 引入:什么是图像? ... -
hslogic_设计匹配QMF和DFT滤波器组的新方法
2020-09-15 18:26:53提出三种不同的滤波器组结构匹配于信号或信号统计信息:2通道均匀滤波器组,M通道并元非均匀滤波器组和M通道改进离散傅里叶变化滤波器组。首先是匹配于信号或信号统计信息的2通道QMF分析滤波器组 ,为了得到这个... -
频率域电磁剖面有限差分法2.5维正演数值模拟
2020-07-13 21:49:52选取适当的kx值,用有限差分法在y-z平面的网格中求解,再通过反傅里叶变换得到空间域中的电磁场。在验证了算法的正确性之后,对不同埋深的直立异常体、倾斜异常体及断陷模型进行了数值模拟,其结果直观地显示了异常体磁... -
基于C#的基础数字图像处理(Visual Studio)
2019-06-26 09:01:51内容有文件(打开、保存、退出)、点处理(彩色转灰阶、图像取反、图像旋转、图像镜像、均衡化、直方图、亮度变化、直方图扩展)、空频域变换(傅里叶变换、傅里叶反变换、离散余弦变换、离散余弦反变换、巴特沃斯低...