应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 帮忙看一下这个程序啊！！
51/1返回列表
查看: 772  |  回复: 4
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

破碎东瓜

金虫 (初入文坛)

[求助] 帮忙看一下这个程序啊！！ 已有1人参与

SUBROUTINE KKFFT(PR,PI,N,K,FR,FI,L,IL)
        DIMENSION PR(N),PI(N),FR(N),FI(N)
        DOUBLE PRECISION PR,PI,FR,FI,P,Q,S,VR,VI,PODDR,PODDI
        DO 20 IT=0,N-1
          M=IT
          IS=0
          DO 10 I=0,K-1
            J=M/2
            IS=2*IS+(M-2*J)
            M=J
10          CONTINUE
          FR(IT+1)=PR(IS+1)
          FI(IT+1)=PI(IS+1)
20        CONTINUE
        PR(1)=1.0
        PI(1)=0.0
        PR(2)=COS(6.283185306/N)
        PI(2)=-SIN(6.283185306/N)
        IF (L.NE.0) PI(2)=-PI(2)
        DO 30 I=3,N
          P=PR(I-1)*PR(2)
          Q=PI(I-1)*PI(2)
          S=(PR(I-1)+PI(I-1))*(PR(2)+PI(2))
          PR(I)=P-Q
          PI(I)=S-P-Q
30        CONTINUE
        DO 40 IT=0,N-2,2
          VR=FR(IT+1)
          VI=FI(IT+1)
          FR(IT+1)=VR+FR(IT+2)
          FI(IT+1)=VI+FI(IT+2)
          FR(IT+2)=VR-FR(IT+2)
          FI(IT+2)=VI-FI(IT+2)
40        CONTINUE
        M=N/2
        NV=2
        DO 70 L0=K-2,0,-1
          M=M/2
          NV=2*NV
          DO 60 IT=0,(M-1)*NV,NV
          DO 60 J=0,(NV/2)-1
            P=PR(M*J+1)*FR(IT+J+1+NV/2)
            Q=PI(M*J+1)*FI(IT+J+1+NV/2)
            S=PR(M*J+1)+PI(M*J+1)
            S=S*(FR(IT+J+1+NV/2)+FI(IT+J+1+NV/2))
            PODDR=P-Q
            PODDI=S-P-Q
            FR(IT+J+1+NV/2)=FR(IT+J+1)-PODDR
            FI(IT+J+1+NV/2)=FI(IT+J+1)-PODDI
            FR(IT+J+1)=FR(IT+J+1)+PODDR
            FI(IT+J+1)=FI(IT+J+1)+PODDI
60          CONTINUE
70        CONTINUE


        IF (L.NE.0) THEN
          DO 80 I=1,N
            FR(I)=FR(I)/N
            FI(I)=FI(I)/N
80          CONTINUE
        END IF
        IF (IL.NE.0) THEN
          DO 90 I=1,N
            PR(I)=SQRT(FR(I)*FR(I)+FI(I)*FI(I))
            PI(I)=ATAN(FI(I)/FR(I))*360.0/6.283185306
90          CONTINUE
        END IF
        RETURN
        END
这个是徐士良的快速傅里叶变换变换算法的源代码，能不能帮忙看一下划线部分到底是什么意思
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
如果你想拥有你从未有过的东西，那么你就要做你从未做过的事情。
1楼 2015-09-15 21:06:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cob

木虫 (正式写手)

【答案】应助回帖

感谢参与，应助指数 +1
do Number xxx
xxxxx
Number continue

循环语句
赞一下(1人) 回复此楼
2楼2015-09-16 20:15:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

破碎东瓜

金虫 (初入文坛)
引用回帖:
2楼: Originally posted by cob at 2015-09-16 20:15:17
do Number xxx
xxxxx
Number continue

循环语句

我想问的是这个部分在整个FFT中起什么作用，那个代码我肯定是懂的，就是不懂为什么要这样编写
赞一下 回复此楼
如果你想拥有你从未有过的东西，那么你就要做你从未做过的事情。
3楼2015-09-17 08:11:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cob

木虫 (正式写手)
引用回帖:
3楼: Originally posted by 破碎东瓜 at 2015-09-17 08:11:06
我想问的是这个部分在整个FFT中起什么作用，那个代码我肯定是懂的，就是不懂为什么要这样编写...

作用应该是找出对应  PR(IS+1) / PI(IS+1) 赋值给 对应 的 FR(IT+1)/FI(IT+1)。
至于为啥这么写，书的公式应该有对应的解释。我没有这本书，所以没办法解释。
赞一下 回复此楼
4楼2015-09-17 19:48:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cob

木虫 (正式写手)
引用回帖:
4楼: Originally posted by cob at 2015-09-17 19:48:03
作用应该是找出对应  PR(IS+1) / PI(IS+1) 赋值给 对应 的 FR(IT+1)/FI(IT+1)。
至于为啥这么写，书的公式应该有对应的解释。我没有这本书，所以没办法解释。...

补充一点，
J=M/2 ， J 是等于 M/2 的整数部分。         
IS=2*IS+(M-2*J)， 所以 M-2*J ，M是奇数等于1，偶数等于0.
赞一下 回复此楼
5楼2015-09-17 20:02:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 破碎东瓜 的主题更新
51/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭