CFUNC.CPP
上传用户:zengbais
上传日期:2022-08-08
资源大小:49k
文件大小:9k
开发平台:

C++ Builder

  1. #include <values.h>
  2. #include <math.h>
  3. #include <stdio.h>
  4. #include "cfunc.h"
  5. COMPLEX calgo::cal(DOUBLE x) { // 基类的基本算法
  6. return yfactor*calculate(xfactor*(x-xshift))+addconst;
  7. }
  8. calgo * calgo::clone() // 克隆自己,必须被继承子类改写
  9. {
  10. return new calgo(*this);
  11. }
  12. calgo * calgo::mul(COMPLEX a) // 乘a
  13. {
  14. calgo * alg;
  15. alg = this;
  16. if(refnum>1) { refnum--;
  17. alg = clone(); // 如引用数大于1,则产生新的对象
  18. }
  19. alg->yfactor *= a;
  20. alg->addconst *= a;
  21. return alg;
  22. }
  23. calgo * calgo::add(COMPLEX a) // 加a
  24. {
  25. calgo * alg;
  26. alg = this;
  27. if(refnum>1) { refnum--;
  28. alg = clone(); // 如引用数大于1,则产生新的对象
  29. }
  30. alg->addconst += a;
  31. return alg;
  32. }
  33. calgo * calgo::neg() // 取负
  34. {
  35. calgo * alg;
  36. alg = this;
  37. if(refnum>1) { refnum--;
  38. alg = clone(); // 如引用数大于1,则产生新的对象
  39. }
  40. alg->addconst = -alg->addconst;
  41. alg->yfactor = -alg->yfactor;
  42. return alg;
  43. }
  44. calgo * calgo::setxfactor(DOUBLE x) // 设置x轴因子
  45. {
  46. calgo * alg;
  47. alg = this;
  48. if(refnum>1) { refnum--;
  49. alg = clone(); // 如引用数大于1,则产生新的对象
  50. }
  51. alg->xfactor = x;
  52. return alg;
  53. }
  54. calgo * calgo::xroom(DOUBLE x) // 将xfactor扩大x倍
  55. {
  56. calgo * c;
  57. if(x != 0)
  58. c = setxshift(xshift/x);
  59. return c->setxfactor(x*xfactor);
  60. }
  61. calgo * calgo::setxshift(DOUBLE x) // 设置xshift的值
  62. {
  63. calgo * alg;
  64. alg = this;
  65. if(refnum>1) { refnum--;
  66. alg = clone(); // 如引用数大于1,则产生新的对象
  67. }
  68. alg->xshift = x;
  69. return alg;
  70. }
  71. calgo * calgo::xshiftas(DOUBLE x) // 从当前开始平移x
  72. {
  73. return setxshift(xshift+x);
  74. }
  75. calgo * calgojoin::clone() // 克隆自己
  76. {
  77. return new calgojoin(*this);
  78. }
  79. COMPLEX calgojoin::calculate(DOUBLE x) // 实施结合算法
  80. {
  81. if(leftalgo==0 || rightalgo==0)
  82. throw TMESSAGE("empty algo pointor!");
  83. COMPLEX a,b;
  84. a = leftalgo->cal(x);
  85. b = rightalgo->cal(x);
  86. if(met == cadd) // 返回各结合运算值
  87. return a+b;
  88. else if(met == csub)
  89. return a-b;
  90. else if(met == cmul)
  91. return a*b;
  92. else if(met == cdiv)
  93. return a/b;
  94. else if(met == cpow)
  95. return pow(a,b);
  96. return 0.0;
  97. }
  98. calgo * calgofun::clone() // 克隆自己
  99. {
  100. return new calgofun(*this);
  101. }
  102. COMPLEX calgofun::calculate(DOUBLE x) // 实施函数算法
  103. {
  104. if(f)
  105. return f(x);
  106. return 0.0;
  107. }
  108. COMPLEX calgoenter::calculate(DOUBLE x) // 实施函数算法
  109. {
  110. return COMPLEX(er->cal(x),ei->cal(x));
  111. }
  112. calgo * calgoenter::clone() // 克隆自己
  113. {
  114. return new calgoenter(*this);
  115. }
  116. cfunc::cfunc()// 缺省构造函数,产生函数f(x)=x;
  117. {
  118. alg = new calgo();
  119. }
  120. cfunc::cfunc(cfunc & fn) // 拷贝构造函数
  121. {
  122. alg = fn.alg;
  123. alg->refnum++;
  124. }
  125. cfunc::cfunc(COMPLEX (*fun)(DOUBLE)) // 函数指针的构造函数
  126. {
  127. alg = new calgofun(fun);
  128. }
  129. cfunc::cfunc(COMPLEX a) // 常函数构造函数
  130. {
  131. alg = new calgo(a);
  132. }
  133. cfunc& cfunc::operator=(cfunc& fn) // 赋值运算符
  134. {
  135. if(this == &fn) return (*this); // 如果等于自己,干脆啥也别做
  136. if(alg) {
  137. alg->refnum--; // 原算法引用数减1
  138. if(!alg->refnum) // 如结果为0则删除原算法
  139. delete alg;
  140. }
  141. alg = fn.alg; // 将fn的算法移过来
  142. if(alg) alg->refnum++; // 引用数增加
  143. return (*this);
  144. }
  145. cfunc& cfunc::operator=(COMPLEX (*fn)(DOUBLE)) // 用函数指针的赋值运算符
  146. {
  147. if(alg) {
  148. alg->refnum--; // 原算法引用数减1
  149. if(!alg->refnum) // 如结果为0则删除原算法
  150. delete alg;
  151. }
  152. alg = new calgofun(fn);
  153. return (*this);
  154. }
  155. cfunc& cfunc::operator=(COMPLEX a) // 常函数的赋值运算符
  156. {
  157. if(alg) {
  158. alg->refnum--; // 原算法引用数减1
  159. if(!alg->refnum) // 如结果为0则删除原算法
  160. delete alg;
  161. }
  162. alg = new calgo(a);
  163. return (*this);
  164. }
  165. cfunc& cfunc::operator+=(cfunc& fn) // 自身加一个函数
  166. {
  167. calgo * a = new calgojoin(alg, fn.alg, cadd);
  168. alg->refnum--; // 因为联合算法对两个算法都加了引用数,因此再减回来
  169. alg = a; // 将新指针赋给alg
  170. return (*this);
  171. }
  172. cfunc& cfunc::operator+=(COMPLEX (*f)(DOUBLE)) // 自身加一个函数指针
  173. {
  174. cfunc fn(f); // 将函数指针包装为函数类
  175. operator+=(fn); // 实施加操作
  176. return (*this);
  177. }
  178. cfunc cfunc::operator+(cfunc& fn) // 相加产生新函数
  179. {
  180. calgo * a = new calgojoin(alg, fn.alg, cadd);
  181. cfunc f(a);
  182. return f;
  183. }
  184. cfunc cfunc::operator+(COMPLEX a) // 与常数相加产生新函数
  185. {
  186. cfunc f(*this);
  187. f += a;
  188. return f;
  189. }
  190. cfunc cfunc::operator+(COMPLEX (*f)(DOUBLE)) // 加一个函数指针产生新函数
  191. {
  192. cfunc ff(*this);
  193. ff += f;
  194. return ff;
  195. }
  196. cfunc& cfunc::neg() // 自身取负
  197. {
  198. alg=alg->neg();
  199. return (*this);
  200. }
  201. cfunc cfunc::operator-() // 产生负函数
  202. {
  203. cfunc f(*this);
  204. f.neg();
  205. return f;
  206. }
  207. cfunc& cfunc::operator-=(cfunc& fn) // 自身减一个函数
  208. {
  209. calgo * a = new calgojoin(alg, fn.alg, csub);
  210. alg->refnum--; // 因为联合算法对两个算法都加了引用数,因此再减回来
  211. alg = a; // 将新指针赋给alg
  212. return (*this);
  213. }
  214. cfunc& cfunc::operator-=(COMPLEX (*f)(DOUBLE)) // 自身减一个函数指针
  215. {
  216. cfunc fn(f); // 将函数指针包装为函数类
  217. operator-=(fn); // 实施减操作
  218. return (*this);
  219. }
  220. cfunc cfunc::operator-(cfunc& fn) // 相减产生新函数
  221. {
  222. calgo * a = new calgojoin(alg, fn.alg, csub);
  223. cfunc f(a);
  224. return f;
  225. }
  226. cfunc cfunc::operator-(COMPLEX a) // 与常数相减产生新函数
  227. {
  228. cfunc f(*this);
  229. f -= a;
  230. return f;
  231. }
  232. cfunc operator-(COMPLEX a, cfunc& f) // 常数减函数
  233. {
  234. return (-f)+a;
  235. }
  236. cfunc cfunc::operator-(COMPLEX (*f)(DOUBLE)) // 减一个函数指针产生新函数
  237. {
  238. cfunc ff(*this);
  239. ff -= f;
  240. return ff;
  241. }
  242. cfunc operator-(COMPLEX (*f)(DOUBLE),cfunc& fn) // 函数指针减函数
  243. {
  244. cfunc ff(f);
  245. ff -= fn;
  246. return ff;
  247. }
  248. cfunc& cfunc::operator*=(cfunc& fn) // 自身乘一个函数
  249. {
  250. calgo * a = new calgojoin(alg, fn.alg, cmul);
  251. alg->refnum--; // 因为联合算法对两个算法都加了引用数,因此再减回来
  252. alg = a; // 将新指针赋给alg
  253. return (*this);
  254. }
  255. cfunc& cfunc::operator*=(COMPLEX (*f)(DOUBLE)) // 自身乘一个函数指针
  256. {
  257. cfunc fn(f); // 将函数指针包装为函数类
  258. operator*=(fn); // 实施乘操作
  259. return (*this);
  260. }
  261. cfunc cfunc::operator*(cfunc& fn) // 相乘产生新函数
  262. {
  263. calgo * a = new calgojoin(alg, fn.alg, cmul);
  264. cfunc f(a);
  265. return f;
  266. }
  267. cfunc cfunc::operator*(COMPLEX a) // 与常数相乘产生新函数
  268. {
  269. cfunc f(*this);
  270. f *= a;
  271. return f;
  272. }
  273. cfunc cfunc::operator*(COMPLEX (*f)(DOUBLE)) // 乘一个函数指针产生新函数
  274. {
  275. cfunc ff(*this);
  276. ff *= f;
  277. return ff;
  278. }
  279. cfunc& cfunc::operator/=(cfunc& fn) // 自身除以一个函数
  280. {
  281. calgo * a = new calgojoin(alg, fn.alg, cdiv);
  282. alg->refnum--; // 因为联合算法对两个算法都加了引用数,因此再减回来
  283. alg = a; // 将新指针赋给alg
  284. return (*this);
  285. }
  286. cfunc& cfunc::operator/=(COMPLEX (*f)(DOUBLE)) // 自身除以一个函数指针
  287. {
  288. cfunc fn(f); // 将函数指针包装为函数类
  289. operator/=(fn); // 实施除法操作
  290. return (*this);
  291. }
  292. cfunc cfunc::operator/(cfunc& fn) // 相除产生新函数
  293. {
  294. calgo * a = new calgojoin(alg, fn.alg, cdiv);
  295. cfunc f(a);
  296. return f;
  297. }
  298. cfunc cfunc::operator/(COMPLEX a) // 与常数相除产生新函数
  299. {
  300. cfunc f(*this);
  301. f /= a;
  302. return f;
  303. }
  304. cfunc operator/(COMPLEX a, cfunc& f) // 常数除以函数
  305. {
  306. cfunc ff(a);
  307. return ff/f;
  308. }
  309. cfunc cfunc::operator/(COMPLEX (*f)(DOUBLE)) // 除以一个函数指针产生新函数
  310. {
  311. cfunc ff(*this);
  312. ff /= f;
  313. return ff;
  314. }
  315. cfunc operator/(COMPLEX (*f)(DOUBLE),cfunc& fn) // 函数指针除以函数
  316. {
  317. cfunc ff(f);
  318. ff /= fn;
  319. return ff;
  320. }
  321. void cfunc::setxfactor(DOUBLE a) // 设置x因子为a
  322. {
  323. alg = alg->setxfactor(a);
  324. }
  325. void cfunc::xroom(DOUBLE a)   // x方向扩大a倍
  326. {
  327. alg = alg->xroom(a);
  328. }
  329. cfunc& cfunc::power(cfunc& f) // 函数的f次乘幂,函数自身改变
  330. {
  331. calgo * a = new calgojoin(alg, f.alg, cpow);
  332. alg->refnum--; // 因为联合算法对两个算法都加了引用数,因此再减回来
  333. alg = a; // 将新指针赋给alg
  334. return (*this);
  335. }
  336. cfunc cfunc::operator^(cfunc& fn) // f次乘幂,产生新函数
  337. {
  338. calgo * a = new calgojoin(alg, fn.alg, cpow);
  339. cfunc f(a);
  340. return f;
  341. }
  342. cfunc& cfunc::power(COMPLEX a) // 函数的a次幂,函数自身改变
  343. {
  344. cfunc f(a);
  345. return power(f);
  346. }
  347. cfunc cfunc::operator^( COMPLEX a)  // 函数的a次幂,产生新函数,原函数不变
  348. {
  349. cfunc f(a);
  350. cfunc ff(*this);
  351. return ff.power(f);
  352. }
  353. COMPLEX calgopoly::calculate(DOUBLE x) // 实施函数算法
  354. {
  355. size_t i;
  356. COMPLEX u;
  357. size_t n=data.rownum-1;
  358. u=data(n);
  359. for (i=n; i>0; i--)
  360. u=u*x+data(i-1);
  361. return u;
  362. }
  363. calgo * calgopoly::clone() // 克隆自己
  364. {
  365. return new calgopoly(*this);
  366. }
  367. cfunc::cfunc(cmatrix& m) // 构造数值相关函数,
  368. // m为nX1矩阵,是n-1阶多项式系数,
  369. // 其中m(0,0)为常数项,m(n-1,0)为n-1次项。
  370. {
  371. alg = new calgopoly(m); // 产生多项式算法
  372. }
  373. void cfunc::setxshift(DOUBLE a) // 设置函数沿x轴平移a
  374. {
  375. alg = alg->setxshift(a);
  376. }
  377. void cfunc::shiftxas(DOUBLE a) // 函数沿x轴右移a
  378. {
  379. alg = alg->xshiftas(a);
  380. }
  381. cfuncenter::cfuncenter(cmatrix& s,DOUBLE t0, DOUBLE dt):
  382. cfunc(new calgoenter(s.real(),s.image(),t0,dt))
  383. {}
  384. cfunc fourier(func& f,DOUBLE tb, DOUBLE te, DOUBLE dt, DOUBLE df)
  385. // 利用fft技术对函数f作傅里叶变换,其中tb为采样窗口的起始点,te为结束点,
  386. // 必须te>tb,dt为采样间隔,df为频率采样间隔,但返回的cfunc是作了插值的
  387. // 插值函数
  388. {
  389. if(tb>=te) throw TMESSAGE("begin point must less than end point");
  390. DOUBLE tspace;
  391. tspace = te-tb;
  392. size_t n;
  393. if(tspace > 1.0/df) df = 1.0/tspace;
  394. n = ceil(1.0/(df*dt));
  395. cmatrix s(n);
  396. s = 0.0;
  397. n = ceil(tspace/dt);
  398. for(size_t i=0; i<n; i++)
  399. s.set(i,f(tb+i*dt));
  400. s.fft();
  401. n = s.rownum;
  402. cmatrix p(n);
  403. df = 1.0/(dt*n);
  404. for(i=0;i<n/2;i++) {
  405. p.set(i,s(i+n/2));
  406. p.set(i+n/2,s(i));
  407. }
  408. for(i=1; i<n/2; i++) {
  409. COMPLEX c;
  410. DOUBLE phi;
  411. phi = 2.0*M_PI*i*df*(-tb);
  412. c = COMPLEX(cos(phi),sin(phi));
  413. p.set(i+n/2,p(i+n/2)*c);
  414. p.set(-i+n/2,p(-i+n/2)*conj(c));
  415. }
  416. p *= dt;
  417. DOUBLE fbegin = -0.5/dt;
  418. return cfuncenter(p,fbegin,df);
  419. }