FFT的应用

FFT的应用

——讲稿

概述

FFT的模板很简单,大家都会背,于是出题的空间就在于建模了。FFT的题目难在建模,往往需要将问题抽象出来,经过一系列转化后得到乘积式的和,再赋予式子各个项的系数一定的意义即可。

一些链接:

FFT快速傅里叶变换

NTT数论变换

MTT:任意模数NTT

分治FFT

多项式求逆

FFT的应用

基本形式

对于类似\(\displaystyle \sum_{i+j=N+k}a_ib_j\)的式子,可以直接通过FFT计算。
其中N是定值,表示元素个数;k是变量,是题目中的系数,a和b是间接已知数组,当a与b的系数和为定值时,可将一个翻转,否则直接计算。

例题:P3723 [AH2017/HNOI2017]礼物

给出循环节长度为\(n\) 的数组\(a, b\),求
\[\min_{c,k}\{\sum_{i=1}^n(c+a_{i+k}-b_i)^2\}\]
其中\(1\leq n \leq 5\times 10^4,1\leq a_i,b_i\leq 100\)

题解:

首先,有一个结论:两个手环增加非负整数亮度,等于其中一个增加一个整数亮度(可以为负)

我们令增加量为\(x\),旋转以后的原数列为\(\{a\}\{b\}\)那么现在的费用就是:
\[\sum_{i=1}^n\left(a_i+x-b_i\right)^2\]
我们把第i项拿出来拆开,得到:
\[\left(a_i+x-b_i\right)^2=a_i^2+b_i^2+x^2+2a_ix-2a_ib_i-2b_ix\]
那么原式变成了
\[\sum_{i=1}^na_i^2+\sum_{i=1}^nb_i^2+nx^2+2x\left(\sum_{i=1}^na_i-\sum_{i=1}^nb_i\right)-2\sum_{i=1}^na_ib_i\]

我们发现,这个式子除了最后一项之外都是确定的,
那么我们只要令最后一项最大,那么就可以得到最小的费用值了现在问题转化为求\[\sum_{i=1}^na_ib_i\]
这个形式非常卷积,我们把数列\(\{a\}\)反过来,变成
\[\sum_{i=1}^na_{n-i+1}b_i\]
这不是一个卷积吗~所以把反过来的数列\(\{a\}\)倍长(为了循环卷积),和数列\(\{b\}\)卷积,得到的项里面的第\(n+1\)\(n\times2\)项的最大值,就是
\[\sum_{i=1}^na_ib_i\]
的最大值然后把前面的不变项加上,就是答案了

代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef long double LD;
const int INF=1e9+7,MAXN=1e6+10;
const LD Pi=3.141592653589793238462643383279;
struct com{
    LD real,imag;
    com(){
        real=imag=0;
    }
    com(LD x,LD y){
        real=x;
        imag=y;
    }
    friend inline com operator+(com x,com y){
        return com(x.real+y.real,x.imag+y.imag);
    }
    friend inline com operator-(com x,com y){
        return com(x.real-y.real,x.imag-y.imag);
    }
    friend inline com operator*(com x,com y){
        return com(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real);
    }
}f[MAXN],g[MAXN];
int N,M,lim=1,L,rev[MAXN];
inline void FFT(com *a,int type){
    for(int i=0;i<lim;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int hf=1;hf<lim;hf<<=1){
        int len=hf<<1;
        com Wn=com(cos(2.0*Pi/len),type*sin(-2.0*Pi/len));
        for(int j=0;j<lim;j+=len){
            com w=com(1,0);
            for(int k=0;k<hf;k++){
                com t1=a[j+k],t2=a[j+k+hf]*w;
                a[j+k]=t1+t2;
                a[j+k+hf]=t1-t2;
                w=w*Wn;
            }
        }
    }
    if(type==-1)
        for(int i=0;i<lim;i++)
            a[i].real=a[i].real/lim+0.5;
}
LL a[MAXN],b[MAXN],sum1,ans,sum2;
inline LL func1(LL x){
    return N*x*x+2*sum1*x;
}
int main(){
    scanf("%d%d",&N,&M);
    for(int i=1;i<=N;i++)
        scanf("%lld",a+i);
    for(int i=1;i<=N;i++)
        scanf("%lld",b+i);
    for(int i=1;i<=N;i++){
        ans+=a[i]*a[i]+b[i]*b[i];
        sum1+=a[i]-b[i];
    }
    ans+=min(func1(floor(-(LD)sum1/N)),func1(ceil(-(LD)sum1/N)));
    /*--------------------convolution---------------------*/
    for(int i=1;i<=N;i++){
        f[i].real=a[N-i+1];
        g[i].real=g[i+N].real=b[i];
    }
    while(lim<=N*3){
        lim<<=1;
        L++;
    }
    for(int i=0;i<lim;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    FFT(f,1);
    FFT(g,1);
    for(int i=0;i<lim;i++)
        f[i]=f[i]*g[i];
    FFT(f,-1);
    for(int i=N;i<=2*N;i++)
        sum2=max(sum2,(LL)f[i].real);
    printf("%lld",ans-2*sum2);
    return 0;
}

直接计算卷积

一些问题经过简单的推导即可推出可以用计算卷积来解决,这类问题多形如对所有不大于\(n\)\(k\)求出某个东西,其中\(k\)的答案为求某个卷积后结果数组的第\(k\)项。

例题:CF993E

给出一个大小为\(n\)的数组\(a\)和一个数\(x\),对于\(0\)\(n\)之间的所有\(k\),求有多少个\(a\)的区间中恰有\(k\)个数小于\(x\)
\(1\leq n\leq 2 \times 10^5, -10^9\leq x,ai \leq 10^9\)

题解

由于\(x\) 固定,直接把\(a\)中小于\(x\) 数的设为\(1\),其余数设为\(0\),只需要求有多少个区间的和为\(k\)。设\(s\)\(a\) 的前缀和,则只需要求有多少个数对\(s_i, s_j\) 满足\(s_j-s_i = k\),或\(s_j - k = s\)i。设\(f_{s_i}\)\(s_i\) 的出现次数,则\(k\) 的答案为
\[\sum_{i=k}^nf_if_{i-k}=\sum_{i=k}^nf_if_{n-i+k}=\sum_{i+j=n+k}f_if_{n-j}\]
\(g_i=f_{n-i}\),求\(f\)\(g\)的卷积即可。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int INF=1e9+7,MAXN=1e6+10/*Min:2^20+10*/;
const long double Pi=3.141592653589793238462643383279;
struct COM{
    long double real,imag;
    COM(){
        real=imag=0;
    }
    COM(long double x,long double y){
        real=x;
        imag=y;
    }
    friend COM operator+(COM x,COM y){
        return COM(x.real+y.real,x.imag+y.imag);
    }
    friend COM operator-(COM x,COM y){
        return COM(x.real-y.real,x.imag-y.imag);
    }
    friend COM operator*(COM x,COM y){
        return COM(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real);
    }
}f[MAXN],g[MAXN];
int N,X,limit=1,lg,rev[MAXN];
inline void FFT(COM *a,int type){
    for(int i=0;i<limit;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int half=1;half<limit;half<<=1){
        int len=half<<1;
        COM Wn=COM(cos(Pi/half),sin(Pi/half)*type);
        for(int j=0;j<limit;j+=len){
            COM w=COM(1,0),x,y;
            for(int k=0;k<half;k++){
                x=a[j+k],y=a[j+k+half]*w;
                a[j+k]=x+y;
                a[j+k+half]=x-y;
                w=w*Wn;
            }
        }
    }
    if(type==-1)
        for(int i=0;i<limit;i++)
            a[i].real/=(double)limit;
}
LL s[MAXN],ans[MAXN];
char ch;
int main(){
    scanf("%d%d",&N,&X);
    for(int i=1;i<=N;i++){
        scanf("%lld",s+i);
        s[i]=(s[i]<X)+s[i-1];
    }
    for(int i=0;i<N;i++){
        f[s[i]].real++;
        g[s[N]-s[i+1]].real++;
    }
    while(limit<=(s[N]<<1)){
        limit<<=1;
        lg++;
    }
    for(int i=0;i<limit;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
    FFT(f,1);
    FFT(g,1);
    for(int i=0;i<limit;i++)
        f[i]=f[i]*g[i];
    FFT(f,-1);
    LL cnt=0;
    for(int i=1;i<=N;i++)
        if(s[i]==s[i-1])
            cnt++;
        else{
            ans[0]+=cnt*(cnt-1)/2+cnt;
            cnt=0;
        }
    ans[0]+=cnt*(cnt-1)/2+cnt;
    for(int i=0;i<s[N];i++)
        ans[s[N]-i]=f[i].real+1e-2;
    for(int i=0;i<=N;i++)
        printf("%lld ",ans[i]);
    return 0;
}

多项式运算

多项式乘法可以用来表示卷积,而借助多项式的性质,可以分析并解决类型更为广泛的问题。其中,最典型的例子是利用生成函数解决组合计数问题,这往往可以简化推导过程,有时还可以借助专用算法优化复杂度。

字符串匹配

代通配符的字符串匹配

\(s,t\)为字符串,其中\(t\)中某些字符是通配符,可以匹配任意字符,求\(s\)\(t\)中的所有匹配的位置。将通配符设为\(0\),其余字符设为非\(0\)的数,则\(s\)\(k\)处匹配当且仅当
\[\sum_{0\leq i < |s|}t_{i+k}(s_i-t_{i+k})^2=0\]
结合计算卷积,很容易算出左式,即可完成匹配

模板:P4173 残缺的字符串

要求匹配两个字符串A,B,两者都有通配符

\[\sum_{0\leq i < |s|}t_{i+k}s_i(s_i-t_{i+k})^2=0\]

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int INF=1e9+7,MAXN=2e6+100;
const double Pi=3.14159265358979323846;
struct com{
    double real,imag;
    com(){
        real=imag=0;
    }
    com(double x,double y){
        real=x;
        imag=y;
    }
    friend inline com operator+(com x,com y){
        return com(x.real+y.real,x.imag+y.imag);
    }
    friend inline com operator-(com x,com y){
        return com(x.real-y.real,x.imag-y.imag);
    }
    friend inline com operator*(com x,com y){
        return com(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real);
    }
    friend inline com operator*(com x,int y){
        return com(x.real*y,x.imag*y);
    }
}f[MAXN],g[MAXN],h[MAXN];
int lim=1,L,rev[MAXN];
inline void FFT(com *a,int type){
    for(int i=0;i<lim;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int hf=1;hf<lim;hf<<=1){
        int len=hf<<1;
        com Wn=com(cos(Pi/hf),type*sin(Pi/hf));
        for(int j=0;j<lim;j+=len){
            com w=com(1,0);
            for(int k=0;k<hf;k++){
                com t1=a[j+k],t2=a[j+k+hf]*w;
                a[j+k]=t1+t2;
                a[j+k+hf]=t1-t2;
                w=w*Wn;
            }
        }
    }
    if(type!=1)
        for(int i=0;i<lim;i++)
            a[i].real=a[i].real/lim;
}
int M,N;
int a[MAXN],b[MAXN],cnt[MAXN];
vector<int> out;
inline void calc(){
    while(lim<=N+M){
        lim<<=1;
        L++;
    }
    for(int i=0;i<lim;i++)
        rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
    
    for(int i=0;i<M;i++)
        f[i].real=a[i]*a[i]*a[i];
    for(int i=0;i<N;i++)
        g[i].real=b[i];
    FFT(f,1); FFT(g,1);
    for(int i=0;i<lim;i++)
        h[i]=f[i]*g[i];
    
    for(int i=0;i<lim;i++)
        f[i]=g[i]=com();
    for(int i=0;i<M;i++)
        f[i].real=a[i]*a[i];
    for(int i=0;i<N;i++)
        g[i].real=b[i]*b[i];
    FFT(f,1); FFT(g,1);
    for(int i=0;i<lim;i++)
        h[i]=h[i]-f[i]*g[i]*2;
    
    for(int i=0;i<lim;i++)
        f[i]=g[i]=com();
    for(int i=0;i<M;i++)
        f[i].real=a[i];
    for(int i=0;i<N;i++)
        g[i].real=b[i]*b[i]*b[i];
    FFT(f,1); FFT(g,1);
    for(int i=0;i<lim;i++)
        h[i]=h[i]+f[i]*g[i];
    FFT(h,-1);
    for(int i=M-1;i<N;i++)
        if(fabs(h[i].real)<0.5)
            out.push_back(i-M+2);
}
char str1[MAXN],str2[MAXN];
int main(){
    scanf("%d%d %s %s",&M,&N,str1,str2);
    for(int i=0;i<M;i++)
        a[M-i-1]=str1[i]=='*'?0:str1[i]-'a'+1;
    for(int i=0;i<N;i++)
        b[i]=str2[i]=='*'?0:str2[i]-'a'+1;
    calc();
    int siz=out.size();
    printf("%d\n",siz);
    for(int i=0;i<siz;i++)
        printf("%d ",out[i]);
    return 0;
}

给出字符串\(s,t\)和非负整数\(d\),求有多少个\(k\),满足对于所有\(s\)的下标\(i\),都存在距离\(k+i\) 不大于\(d\)\(j\),使得\(s_i = t_j\)\(1\leq |s|,|t|, k\leq 2\times 10^5\),字符集大小为\(4\)

题解

考虑分别处理每种字符。处理字符\(c\) 时,将\(t\) 中所有和字符为\(c\)的位置距离不超过\(d\) 的位置设为\(1\),其余位置设为\(0\);将\(s\) 中所有字符为\(c\) 的位置设为\(1\),其余位置设为\(0\)。那么在\(k\) 处字符\(c\)的匹配位置个数为
\[\sum_{0\leq i\leq|s|}s_it_{i+k}=\sum_{i+j=|s|+k}{t_is_j}\]

可以通过卷积计算。之后只要判断在每个\(k\) 处的每种字符的匹配位置个数之和是否等于\(|s|\) 即可。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef long double LD;
const int INF=1e9+7,MAXN=1e6+10;
const LD Pi=3.141592653589793238462643383279;
struct com{
    LD real,imag;
    com(){
        real=imag=0;
    }
    com(LD x,LD y){
        real=x;
        imag=y;
    }
    friend inline com operator+(com x,com y){
        return com(x.real+y.real,x.imag+y.imag);
    }
    friend inline com operator-(com x,com y){
        return com(x.real-y.real,x.imag-y.imag);
    }
    friend inline com operator*(com x,com y){
        return com(x.real*y.real-x.imag*y.imag,x.real*y.imag+x.imag*y.real);
    }
}f[MAXN],g[MAXN];
int lim=1,L,rev[MAXN];
inline void FFT(com *a,int type){
    for(int i=0;i<lim;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int hf=1;hf<lim;hf<<=1){
        int len=hf<<1;
        com Wn=com(cos(2.0*Pi/len),type*sin(-2.0*Pi/len));
        for(int j=0;j<lim;j+=len){
            com w=com(1,0);
            for(int k=0;k<hf;k++){
                com t1=a[j+k],t2=a[j+k+hf]*w;
                a[j+k]=t1+t2;
                a[j+k+hf]=t1-t2;
                w=w*Wn;
            }
        }
    }
    if(type==-1)
        for(int i=0;i<lim;i++)
            a[i].real=a[i].real/lim+0.5;
}
int N,M,K,ans,cnt[MAXN];
char s[MAXN],t[MAXN];
inline void calc(char ch){
    memset(f,0,sizeof(f));
    memset(g,0,sizeof(g));
    for(int i=0,l=-INF;i<N;i++){
        if(s[i]==ch)
            l=i;
        if(i-l<=K)
            f[i].real=1;
    }
    for(int i=N-1,l=INF;i>=0;i--){
        if(s[i]==ch)
            l=i;
        if(l-i<=K)
            f[i].real=1;
    }
    for(int i=0;i<M;i++)
        g[i].real=t[M-i-1]==ch;
    FFT(f,1);
    FFT(g,1);
    for(int i=0;i<lim;i++)
        f[i]=f[i]*g[i];
    FFT(f,-1);
    for(int i=0;i<N;i++)
        cnt[i]+=floor(f[i+M-1].real);
}
int main(){
    scanf("%d%d%d %s %s",&N,&M,&K,s,t);
    while(lim<=N+M-2){
        lim<<=1;
        L++;
    }
    for(int i=0;i<lim;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for(int i=0;i<4;i++)
        calc("AGCT"[i]);
    for(int i=0;i<N;i++)
        ans+=cnt[i]>=M;
    printf("%d",ans);
    return 0;
}

转载于:https://www.cnblogs.com/guoshaoyang/p/11296027.html

  • 1
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
HAL库是针对STM32系列微控制器的硬件抽象层库。它提供了一组函数和驱动程序,用于简化对STM32微控制器的配置和操作。ADC(模数转换器)是一种用于将模拟信号转换为数字信号的设备。FFT(快速傅里叶变换)是一种用于将时域信号转换为频域信号的算法。 在使用HAL库进行ADC和FFT应用时,可以按照以下步骤进行操作: 1. 配置ADC模块:使用HAL库提供的函数设置ADC的采样时间和分辨率\[1\]。 2. 启动ADC转换:使用HAL库提供的函数启动ADC转换,并使用DMA(直接内存访问)方式将采样数据存储到指定的缓冲区中\[2\]。 3. 等待转换完成:使用一个标志位(例如AdcConvEnd)来检测ADC转换是否完成\[3\]。 4. 进行FFT变换:使用FFT算法对采样数据进行频域转换,将时域信号转换为频域信号。 5. 分析和处理频域数据:根据应用需求,对FFT得到的频域数据进行分析和处理,例如频谱分析、滤波等。 需要注意的是,具体的代码实现可能会因为使用的STM32系列微控制器型号和HAL库版本而有所不同。因此,在实际应用中,需要参考相关的文档和示例代码来进行具体的配置和操作。 希望以上信息对您有所帮助! #### 引用[.reference_title] - *1* [STM32 CubeMX配置ADC+DMA进行FFT(1)](https://blog.csdn.net/ANNOLRA/article/details/119206287)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v91^insertT0,239^v3^insert_chatgpt"}} ] [.reference_item] - *2* *3* [STM32HAL ADC+TIM+DMA采集交流信号 基于cubemx](https://blog.csdn.net/qq_34022877/article/details/121941236)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v91^insertT0,239^v3^insert_chatgpt"}} ] [.reference_item] [ .reference_list ]

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值