houdini fft模型

一路向难走

已于 2024-01-21 14:33:22 修改

阅读量402

点赞数 6

文章标签： houdini

于 2023-12-29 14:56:29 首次发布

本文链接：https://blog.csdn.net/qq_39239990/article/details/135290459

版权

文章讨论了在Houdini中使用复杂结构进行FFT计算的问题，指出Python的numpy.fft封装导致无法查看底层代码。作者提到使用vector代替complex结构以及尝试了Pffft库和FFTPACK源码，以解决精度和准确性问题。

摘要由CSDN通过智能技术生成

按照博友给的代码，进行可视化，发现绘出的图形没有错误，但计算的结果16，32，64.。。居然是错的；哪位知道正确的，敬请指正，或是后续啥时候复习完后再改正吧

houdini 上无法定义结构，可以使用vector 代替complex结构进行计算

使用网上的fft代码，可以计算，但无法产生理想效果；精度与准确度无法实现

python.numpy.fft可以实现相应的准确度，但是他对fft的底层进行了封装，无法查看

fft的源始代码：FFTPACK

pffft的代码 Bitbucket 此处可获取fftpack.c .h文件（教程：Pffft库的学习-CSDN博客）

主要函数代码



i@c_lay = detail(1,'iteration',0);
//printf("current _lay: %d ",@c_lay);

/*
float wn_r[@len];
float wn_i[@len];

//init wn
for(int k=0;k<@len;k++)
{
    wn_r[k] = cos(2*3.14/@len*k);
    wn_i[k] = -1*sin(2*3.14/@len*k);
    
}
*/



int len = @numpt;
i@t_lay = @lays;  //log2 8
//printf("%d",@numpt);



i@dist = shl(1,@c_lay);        //计算每一级的间隔距离
 
for( int j=0;j<@dist;j++ )      
{

    for(int k=j;k<len;k+=(shl(@dist,1)))   //k计算的序号1   
    {
        int q = k+@dist; //q计算的序号2
        int z = shl(k, (@t_lay - @c_lay -1));    //确定旋转因子的指数
        
        //w = rotation_factor(len,1,z);

        f@w_real = cos(2* PI * z / len);
        f@w_img = -sin(2* PI * z / len);    
        
        f@q_real = point(0,'real',q);
        f@q_img = point(0,'img',q);
        
        f@k_real = point(0,'real',k);
        f@k_img = point(0,'img',k);
     //   printf("k_real for %d ",@c_lay);
     //   printf("%f  ",@k_real);
                        
        float tmp_real =@k_real;
        float tmp_img = @k_img;
        
        //complex_multi x
        f@wxq_real = @w_real*@q_real -@w_img*@q_img;
        f@wxq_img = @w_real*@q_img + @w_img*@q_real;
        
        
        f@wxq_real = @wxq_real * @q_real - @wxq_img * @q_img;
        f@wxq_img =  @wxq_real * @q_img + @wxq_img * @q_real;
      
        
        
        
        //complex_add  x[k]
        tmp_real = @k_real + @wxq_real;
        tmp_img  = @k_img +  @wxq_img;
        
        setpointattrib(0,'real',k,tmp_real);
        setpointattrib(0,'img',k,tmp_img);
        
        //printf("%f",tmp_img);
        
        //complex sub  x[q]
        tmp_real =@k_real-@wxq_real;
        tmp_img = @k_img -@wxq_img;
        
        setpointattrib(0,'real',q,tmp_real);
        setpointattrib(0,'img',q,tmp_img);
        
        //printf("%f",tmp_img);
         //add prim
         
         int plen = detail(0,"len",0);
              v@pk1 = point(0,'P',plen+k);
              v@pq1 = point(0,'P',plen+q);
              
              v@pk2 = point(0,'P',k);
              v@pq2 = point(0,'P',q);  
              
          //    i@p1 = addpoint(0,@pk1);
          //    i@p2 = addpoint(0,@pq1);
              
              addprim(0,'polyline',k,plen+q);
              addprim(0,'polyline',q,plen+k);
        
       

    }
}

每次循环中进行两次计算（一次算两）

8点 = log2(8) = 3 分层；32点 5层，64点 6层

每层点的跨度 dist << n(n代表层）

python 直接实现fft计算

import numpy as np
node = hou.pwd()
geo = node.geometry()

# Add code to modify contents of geo.
# Use drop down menu to select examples.
#n = len(geo.points())
#t = geo.points()[0]
#t.setPosition((0,1,0))

a =[]
for i in range(len(geo.points())):
    print(i)
    p = (geo.points()[i].position()[1])
    a.append(p)

ft = np.fft.fft(a)
print(ft[0])  

#aa = geo.addAttrib(hou.attribType.Point,'mreal',0.0)
#cd = geo.addAttrib(hou.attribType.Point,'Cd',(0,1,0))
for i in range(len(geo.points())):
    print(i)
    #p = geo.points()[i].setPosition((0,1,1))
    #geo.points()[i].setAttribValue(aa,0.5)

    geo.points()[i].setAttribValue('real',ft[i].real)
    geo.points()[i].setAttribValue('img',ft[i].imag)

houdini 中使用vector 代替 complex结构进行计算

///-------------------------------------------------------------------------------------------------------------------------

SWAP 代码


int i=0,j=0,k=0;
float t;
int size_x = @numpt;
for(i=0;i<size_x;i++)
{
    k=i;j=0;
    t=(log(size_x)/log(2));
    while( (t--)>0 )    //利用按位与以及循环实现码位颠倒
    {
        j = shl(j,1); //j=j<<1;
        j|=(k & 1);
        k = shr(k,1); //k=k>>1;
    }
    if(j>i)    //将x(n)的码位互换
    {
       // f@ni = point(0,'real',i);
        //f@nj = point(0,'real',j);
        
        //setpointattrib(0,'real',i,@nj);
        //setpointattrib(0,'real',j,@ni);
        v@xi = point(0,'cmplex',i);
        v@xj = point(0,'cmplex',j);
        
        setpointattrib(0,'cmplex',i,@xj);
        setpointattrib(0,'cmplex',j,@xi);
        
        
        //temp=x[i];
        //x[i]=x[j];
        //x[j]=temp;
    }
}

//使用vector 代替 complex 代码 

//int size_x = @numpt;
int size_x = 47;

void xadd(vector a,b,c)
{
        c.x = a.x+b.x;   //x->real ; y->img
        c.y = a.y+b.y;
}

void xmul(vector a,b,c)
{
        c.x = a.x*b.x-a.y*b.y;
        c.y = a.x*b.y+a.y*b.x;
}

void xsub(vector a,b,c)
{
        c.x = a.x-b.x;
        c.y = a.y-b.y;
}

vector W2[size_x];
v[]@W;
for(int i=0;i<size_x;i++)
{
    W2[i].x=cos(2*PI/size_x*i);   //用欧拉公式计算旋转因子
    W2[i].y=-1*sin(2*PI/size_x*i);
}
printf("------w------");
printf("%f",W2);
printf("-------");

        
        
int i=0,j=0,k=0,l=0;
//complex up,down,product;
vector up,down,product;

//change();  //调用变址函数
for(i=0;i< log(size_x)/log(2) ;i++)        /*一级蝶形运算 stage */
{   
        l = shl(1,i);  //l=1<<i;
        for(j=0;j<size_x;j+= 2*l )     /*一组蝶形运算 group,每个group的蝶形因子乘数不同*/
        {            
                for(k=0;k<l;k++)        /*一个蝶形运算 每个group内的蝶形运算*/
                {       
                //f@real = point(0,'real',j+k+l);
                //f@img = point(0,'real',j+k+l);
                v@Xplx = point(0,'cmplex',j+k+l);
                xmul(@Xplx,W2[size_x*k/2/l],product);                
                v@Xplx = point(0,'cmplex',j+k);
                xadd(@Xplx,product,up);
                xsub(@Xplx,product,down);
                
       
                
                setpointattrib(0,'cmplex',j+k,up);
                setpointattrib(0,'cmplex',j+k+l,down);
                
                        //xmul(x[j+k+l],W2[size_x*k/2/l],product);
                        //xadd(x[j+k],product,up);
                        //xsub(x[j+k],product,down);
                        //x[j+k]=up;
                        //x[j+k+l]=down;
                }
        }
}

使用visual studio 2017 中使用 fftpack

Bitbucket 此处获取fftpack.c .h 源码

1.将代码复制到项目文件夹下（****）

测试代码

// ConsoleApplication1.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//

#include "pch.h"
#include <iostream>
#include "fftpack.h"



// macos 转 高级保存选项
// 属性 c++/c 预编辑头文件 否 
//文件复制到 项目 目录下
int main()
{
	int N = 16;

	int Nfloat = 2 * N;  

	int Nbytes = Nfloat * sizeof(fftpack_real);

	fftpack_real *wrk = (fftpack_real*)malloc(2 * Nbytes + 15 * sizeof(fftpack_real));

	fftpack_real *refin = (fftpack_real*)malloc(Nbytes);





	printf("==================\ncomplex orgin data:\n");

	int i = 0;
	for (int k = 0; k < Nfloat; k += 2) {
		
		//refin[k] = k * 0.5f;
		//**********************
		refin[k] = sin(6.28/(N-1)*i);
		printf("%f  ",  6.28/(N-1)*i);
		i = i + 1;

		refin[k + 1] = 0 * k*0.5f;

		printf("%d %f %f\n", k / 2, refin[k], refin[k + 1]);

	}

	cffti(N, wrk);
	printf("------------wrk----------\n");
	
	for (int i = 0; i < Nfloat+15; i++)
	{
		printf("%d  %f \n",i, wrk[i]);
	}

	cfftf(N, refin, wrk);

	for (int i = 0; i < Nfloat + 15; i++)
	{
		printf("%d  %f \n", i, refin[i]);
	}


}

// 运行程序: Ctrl + F5 或调试 >“开始执行(不调试)”菜单
// 调试程序: F5 或调试 >“开始调试”菜单

// 入门提示: 
//   1. 使用解决方案资源管理器窗口添加/管理文件
//   2. 使用团队资源管理器窗口连接到源代码管理
//   3. 使用输出窗口查看生成输出和其他消息
//   4. 使用错误列表窗口查看错误
//   5. 转到“项目”>“添加新项”以创建新的代码文件，或转到“项目”>“添加现有项”以将现有代码文件添加到项目
//   6. 将来，若要再次打开此项目，请转到“文件”>“打开”>“项目”并选择 .sln 文件

结果与python 一致