FFT专题

最新推荐文章于 2024-03-19 14:11:16 发布
最新推荐文章于 2024-03-19 14:11:16 发布
阅读量97 收藏
点赞数 1
分类专栏： FFT 文章标签： c语言 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_52383696/article/details/120813165
版权

FFT

最近在学一个叫做FFT的东西(多项式我来啦

黑暗爆炸 - 4503
模板字符串匹配，卷一下就好啦

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int maxn=3e5+10; 
const double pi=acos(-1.0);

ll n,m;
char ch[maxn];
struct comp
{
    double x,y;
    comp(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn],a2[maxn],b2[maxn],b3;
comp operator + (comp a,comp b){ return comp(a.x+b.x,a.y+b.y);}
comp operator - (comp a,comp b){ return comp(a.x-b.x,a.y-b.y);}
comp operator * (comp a,comp b){ return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

int l=0,r[maxn];
int limit=1;
void FFT(comp a[],int type){
    for(int i=0;i<limit;i++) 
        if(i<r[i])swap(a[i],a[r[i]]);  //求出要迭代的序列 
    for(int mid=1;mid<limit;mid<<=1)  //待合并区间的中点
    {
        comp wn(cos(pi/mid),type*sin(pi/mid));  //单位根 
        for(int R=mid<<1,j=0;j<limit;j+=R)  //R是区间的右端点，j表示前已经到哪个位置了 
        {
            comp w(1,0); //幂 
            for(int k=0;k<mid;k++,w=w*wn)  //枚举左半部分 
            {
                comp x=a[j+k],y=w*a[j+mid+k]; //蝴蝶变化 
                a[j+k]=x+y;
                a[j+mid+k]=x-y;
            }
        }
    }
}
vector<int>v;
int main()
{
	scanf("%s",ch);
	n=strlen(ch);
	for(int i=0;i<n;i++)
		if(ch[i]=='?')a[i].x=0;	 
		else a[i].x=ch[i]-'a'+1;
	scanf("%s",ch);
	m=strlen(ch);
	for(int i=0;i<m;i++)
		if(ch[i]=='?')b[i].x=0;	 
		else b[i].x=ch[i]-'a'+1;
	reverse(b,b+m);
	while(limit<=n+m)limit<<=1,l++;
	for(int i=0;i<limit;i++)
		a2[i]=a[i]*a[i],
		b2[i]=(comp){2,0}*b[i]*b[i],b3.x+=b[i].x*b[i].x*b[i].x;
    for(int i=0;i<limit;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a,1);FFT(b,1);FFT(a2,1);FFT(b2,1);
    for(int i=0;i<limit;i++){
    	a2[i]=a2[i]*b[i]-a[i]*b2[i];
	}
	FFT(a2,-1);
    for(int i=m-1;i<n;i++){
    	if((int)(a2[i].x/limit+b3.x+0.5)==0){
    		v.push_back(i-m+1);
		}
	}
	printf("%d\n",v.size());
	for(int i=0;i<v.size();i++){
		printf("%d\n",v[i]);
	}
	return 0;
}

黑暗爆炸 - 4259
同上题，模板

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int maxn=7e5+10; 
const double pi=acos(-1.0);

ll n,m;
char ch[maxn];
struct comp
{
    double x,y;
    comp(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn],a2[maxn],b2[maxn],b3[maxn],a3[maxn];
comp operator + (comp a,comp b){ return comp(a.x+b.x,a.y+b.y);}
comp operator - (comp a,comp b){ return comp(a.x-b.x,a.y-b.y);}
comp operator * (comp a,comp b){ return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

int l=0,r[maxn];
int limit=1;
void FFT(comp a[],int type){
    for(int i=0;i<limit;i++) 
        if(i<r[i])swap(a[i],a[r[i]]);  //求出要迭代的序列 
    for(int mid=1;mid<limit;mid<<=1)  //待合并区间的中点
    {
        comp wn(cos(pi/mid),type*sin(pi/mid));  //单位根 
        for(int R=mid<<1,j=0;j<limit;j+=R)  //R是区间的右端点，j表示前已经到哪个位置了 
        {
            comp w(1,0); //幂 
            for(int k=0;k<mid;k++,w=w*wn)  //枚举左半部分 
            {
                comp x=a[j+k],y=w*a[j+mid+k]; //蝴蝶变化 
                a[j+k]=x+y;
                a[j+mid+k]=x-y;
            }
        }
    }
}
vector<int>v;
int main()
{
	scanf("%d%d",&n,&m);
	scanf("%s",ch);
	for(int i=0;i<n;i++)
		if(ch[i]=='*')a[i].x=0;	 
		else a[i].x=ch[i]-'a'+1;
	scanf("%s",ch);
	m=strlen(ch);
	for(int i=0;i<m;i++)
		if(ch[i]=='*')b[i].x=0;	 
		else b[i].x=ch[i]-'a'+1;
	reverse(a,a+n);
	while(limit<=n+m)limit<<=1,l++;
	for(int i=0;i<limit;i++)
		a2[i]=a[i]*a[i],
		b2[i]=(comp){2,0}*b[i]*b[i],
		a3[i]=a[i]*a[i]*a[i],
		b3[i]=b[i]*b[i]*b[i];
    for(int i=0;i<limit;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a3,1);FFT(b3,1);FFT(a2,1);FFT(b2,1);
    FFT(a,1);FFT(b,1);
    for(int i=0;i<limit;i++){
    	a2[i]=a3[i]*b[i]-a2[i]*b2[i]+b3[i]*a[i];
	}
	FFT(a2,-1);
    for(int i=n-1;i<m;i++){
    	if((int)(a2[i].x/limit+0.5)==0){
    		v.push_back(i-n+1);
		}
	}
	printf("%d\n",v.size());
	for(int i=0;i<v.size();i++){
		printf("%d ",v[i]+1);
	}
	return 0;
}

CodeForces - 528D
模糊字符串匹配，我太菜了，看了题解。
原来的想法，是预处理出A串中能匹配的字符，然后只要和B串卷积连乘就行了，比如式子是这样的：

但是这个式子展开也太离谱了！
题解的思路：
因为字符集只有’A,C,G,T’，我们把每个字符分开做贡献，如果看从一个位置开始能匹配几个给定字符，然后再把所有的字符累加即可，如果一项的系数是B串的长度，即意味着匹配成功。

对于一个字符c，
如果A串的第i个位置可以放c，那么 f i = 1 fi=1 fi=1，否则 f i = 0 fi=0 fi=0。　
如果B串的第i个位置可以放c，那么 g i = 1 gi=1 gi=1，否则 g i = 0 gi=0 gi=0。
然后卷积：
h i = ∑ f j hi=∑fj hi=∑fj g i − j gi−j gi−j
然后把4个字符的解累加起来。
r e s i = ∑ h [ c ] resi=∑h[c] resi=∑h[c]
显然，如果 r e s [ i ] = ∣ B ∣ res [i]=|B| res[i]=∣B∣，那么从A串的第 i − ∣ B ∣ + 1 i−|B|+1 i−∣B∣+1位开始可以匹配。
然后时间复杂度还是O( ( ∣ A ∣ + ∣ B ∣ ) l o g ( ∣ A ∣ + ∣ B ∣ ) (|A|+|B|)log(|A|+|B|) (∣A∣+∣B∣)log(∣A∣+∣B∣))。常数小了很多，也不用自己展开了。

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int maxn=7e5+10; 
const double pi=acos(-1.0);

ll n,m,k;
struct comp
{
    double x,y;
    comp(double xx=0,double yy=0){x=xx,y=yy;}
}x[maxn],y[maxn],z[maxn];
comp operator + (comp a,comp b){ return comp(a.x+b.x,a.y+b.y);}
comp operator - (comp a,comp b){ return comp(a.x-b.x,a.y-b.y);}
comp operator * (comp a,comp b){ return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

int l=0,r[maxn];
int limit=1;
void FFT(comp a[],int type){
    for(int i=0;i<limit;i++) 
        if(i<r[i])swap(a[i],a[r[i]]);  //求出要迭代的序列 
    for(int mid=1;mid<limit;mid<<=1)  //待合并区间的中点
    {
        comp wn(cos(pi/mid),type*sin(pi/mid));  //单位根 
        for(int R=mid<<1,j=0;j<limit;j+=R)  //R是区间的右端点，j表示前已经到哪个位置了 
        {
            comp w(1,0); //幂 
            for(int k=0;k<mid;k++,w=w*wn)  //枚举左半部分 
            {
                comp x=a[j+k],y=w*a[j+mid+k]; //蝴蝶变化 
                a[j+k]=x+y;
                a[j+mid+k]=x-y;
            }
        }
    }
}
int mp[300],a[maxn],b[maxn];
int res[maxn];
char ch[maxn];
int main()
{
	mp['A']=0;mp['C']=1;mp['G']=2;mp['T']=3;
	scanf("%d%d%d",&n,&m,&k);
	scanf("%s",ch);
	for(int i=0;i<n;i++){
		a[i]=mp[ch[i]];
	}
	scanf("%s",ch);
	for(int i=0;i<m;i++){
		b[i]=mp[ch[i]];
	}
	reverse(b,b+m);
	while(limit<=n+m)limit<<=1,l++;
	for(int i=0;i<limit;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    memset(res,0,sizeof res);
    for(int i=0;i<4;i++){
    	for(int j=0;j<limit;j++){
    		x[j].x=y[j].x=0;
    		x[j].y=y[j].y=0;
		}
		int cnt=0;
		for(int j=0;j<k;j++)
			if(a[j]==i)cnt++;
		for(int j=0;j<n;j++){
			if(j+k<n)
				if(a[j+k]==i)cnt++;
			if(j-k-1>=0)
				if(a[j-k-1]==i)cnt--;
			if(cnt>0)x[j].x=1;
		}
		for(int j=0;j<m;j++){
			if(b[j]==i)y[j].x=1;
		}
		FFT(x,1);FFT(y,1);
		for(int j=0;j<limit;j++){
			z[j]=x[j]*y[j];
		}
		FFT(z,-1);
		for(int j=0;j<n;j++){
			res[j]+=(int(z[j].x/limit+0.5));
		}
	}
	int ans=0;
	for(int i=m-1;i<n;i++)
		if(res[i]==m)ans++;
	printf("%d\n",ans);
	return 0;
}
dplovetree
关注
  • 1
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
fft函数matlab_fft
weixin_39603823的博客
12-04 847
DTFT是离散时间非周期信号的傅里叶变换，由于其频域曲线是连续周期的，计算机无法处理，所以用冲击函数对其频谱在频域进行采样，得到时域和频域都离散的信号，这类过程叫做DFT。写出DFT的公式我们可以看到，完整的求出信号的DFT需要计算N*N次乘法和N*(N-1)次加减法，所以有大牛又提出了FFT算法简化计算过程，加快计算速度。 matlab提供的fft算法有几...
FFT专题训练
秋叶红于二月花
05-14 435
训练的目的是学会如何使用fft解决问题,参考资料 郭晓旭《Fourier transform》 luogu P3083 【模板】多项式乘法（FFT） 直接读入然后求卷积即可,注意精度问题,最开始PI设为3.1415926,WA了三个点,改成4*atan(1.0)就AC了. SPOJ 8372 TSUM - Triple Sums 容斥,FFT 给n(40000)个数(±20000)...
信号处理之快速傅里叶变换(FFT)-通俗易懂
最新发布
红叶幽香的博客
03-19 2万+
历史溯源 欧拉公式 傅里叶级数(FS) 傅里叶变换(FT) 离散傅里叶级数(DFS) 离散时间傅里叶变换(DTFT) 离散傅里叶变换(DFT) 快速傅里叶变换(FFT) MATLAB中常用的FFT工具 FFT中常见的问题(频率混叠，频谱泄露，栅栏效应，旁瓣效应)
FFT/NTT/FWT专题
Qianyri的博客
09-08 557
FFT可快速计算卷积，卷积为两个生成函数的乘积，生成函数为 FFT 卷积 生成函数 的关系 多项式 之 快速傅里叶变换(FFT)/数论变换(NTT)/例题与常用套路【入门】 HDU1402 A * B Problem Plus 大数乘法 生成函数:,设两个大数位：、，L1/L2为两个大数的长度， 利用FFT求两个函数的卷积，将所得改为十进制整数即为所求答案 #include&l...
FFT（多项式）专题水分小记（FFT注意事项）
一位蒟蒻的小博客
05-04 422
DescribePty继续着他的疯狂奔跑，终于渐渐体力不支，在一个应该拐弯的地方没有刹住车，掉入了深深的沼泽中，“啊~~~~~·”pty惊恐的大叫，突然从梦中惊醒了。哪里还有什么奇怪的金字塔，沼泽地，大树。。。只是一个梦而已呀。看了看自己熟悉的房间，pty定了定神。好不容易恢复了过来，pty突然想到还有集训队的互测题没有出！！，如果没有出完的话，后果= =。。啧啧。。pty宁愿再回到金字塔去。于是p
Matlab fftshift 详解
热门推荐
魔法森林博客
05-08 12万+
说明：本文为转载加上自己的编辑，转载请注明出处。 Matlab fftshift 详解 一.实信号情况 因为实信号以fs为采样速率的信号在 fs/2处混叠，所以实信号fft的结果中前半部分对应[0, fs/2],后半部分对应[ -fs/2, 0] 1）实信号fft的结果前半部分对应[0, fs/2]是正频率的结果,后半部分对应[ -fs/2, 0]是负频率的结果。大于fs/2的部分的频谱实...
FFT专题介绍（中文）
02-23
专讲Fourier变换的一部分——FFT的原理、应用、计算机实现，大约25页
专题资料（2021-2022年收藏）基于MATLAB的谐波分析FFT.doc
10-06
最新完美版教育资料
ECG.rar_ecg analysis_心电 向量机_心电分析
09-20
计算机医学应用已成为现代医学的一个重要组成部分和专门学科。计算机化的心电图（computerizedECG）分析是...近年来，欧洲心脏病学大会，国际心电图和向量图学术大会也把心电图的计算机分析列为大会交流的专题之一。
matlab反应扩散代码-SP_Lib:信号处理方法及算法库
05-19
fft专题 滤波 无限长单位脉冲响应（Infinite impulse response,IIR）滤波器 有限长单位脉冲响应（Finite impulse response,FIR）滤波器 由于后续算法变分模态分解中使用到。 扩展的关于Vold-Kalman Order Filtering ...
LTE专题研究之物理层OFDM/SC-FDMA原理
06-17
正交频分复用(OFDM: Orthogonal Frequency Division Multiplexing)是一种多载波调制技术，早在20世纪60年代就已经提出了OFDM的概念，不过由于实现复杂度高，大家并不怎么关注，之后随着DFT(离散傅立叶变化)、FFT...
fft函数C语言实现
12-05
快速傅里叶变换（fft）函数用C语言实现的源代码
C语言FFT函数库及安装教程(Windows和Linux)和使用教程
01-12
压缩包包含：1.C语言FFT函数库FFTW，FFTW 是由麻省理工学院计算机科学实验室超级计算技术组开发的一套离散傅立叶变换(DFT)的计算库，开源、高效和标准 C 语言编写的代码使其得到了非常广泛的应用， Intel 的数学库和 Scilib(类似于 Matlab 的科学计算软件)都使用 FFTW 做 FFT 计算。FFTW 是计算离散 Fourier 变换(DFT)的快速 C 程序的一个完整集合。 2.Windows和Linux下的安装教程和使用教程
fft函数实用实例
06-18
使用fft做常规信号的傅立叶变换，并画图设置坐标轴，横坐标为实际频率，纵坐标取分贝
基于matlab的电力系统谐波仿真
08-20
基于matlab的电力系统谐波仿真，基于matlab的电力系统谐波仿真
Matlab中fft函数的用法及关键问题详解
那年的河川
08-30 11万+
FFT是Fast Fourier Transform(快速傅里叶变换)的简称，这种算法可以减少计算DFT（离散傅里叶变换，关于此更详细的说明见后文）的时间，大大提高了运算效率，并曾经一度被认为是信号分析技术划时代的进步，其重要性由此可见一斑。闲话少叙，言归正传。 基于FFT在信号分析中的重要性，其必然会成为MATLAB的座上宾。FFT算法在MATLAB中实现的函数是Y=fft(x,N)。刚接触频
fft()函数详解
pldcat的博客
01-01 1万+
fft()函数简单到发指，一般使用时就两个参数fft(nparray,n),n还可以缺省。上代码： import numpy as np from scipy.fftpack import fft,ifft fft_y=fft(y) print(fft_y) 执行结果： [180444.84 -0.j -1764.15187386-6325.24578909j 2924.15553426-4550.04232005j ... 733.09073339+342
MATLAB实现的FFT函数
幻雪神界
04-14 3090
写于09年08月29日，函数功能，用MATLAB实现的FFT函数。   完整程序 function xn=myfft(xn,N,M) mybitrevorder(xn,N); for m=1:M B=2^(m-1); for r=0:B-1 P=2^(M-m)*r; for k=r:2^m :N-1 f= xn(k+1); xn(k+
FFT函数的实现(附Matlab代码)
flypassion的博客
06-21 1万+
利用fft函数，实现对100Hz和200Hz信号叠加后信号的频谱分析，分别绘制出了时域和频率的波形。 % fft module % fft function used % Date:2020.06.21 % Author:flypassion clear all close all %%------------------ parameter ------------------%% fs = 2e3; % sampling frequence f1 = 1e2; % signal 1 frequenc
pytorch FFT
09-19
在1.7.1版本之后的PyTorch中，你可以使用`torch.fft.fft`函数进行一维FFT变换，使用`torch.fft.fftn`函数进行n维FFT变换。相应地，你可以使用`torch.fft.ifft`函数进行一维逆FFT变换，使用`torch.fft.ifftn`函数进行...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值