同学你好，我在百度上看到你有球形粒子mie散射相函数matlab源程序，可以给我发一下不，感激不尽！！！

1.在[4]下载message.dat数据文件。文件包含一个带有频率为400HZ的正弦厅扮脊噪音的声波，样本的频率为22000HZ。请按照如下的步骤创建

一个最小二乘的滤波器回复初始的声波。

（1）利用如下的matlab代码来观察文件的频谱：[freq_resp,freq_index]=freqz(message,1,50000,22000)plot(freq_index,abs(freq_resp))

（2）你会看到在频域下文件经过傅里叶变换后的最大值点。这时，你能看出来包含在文件中的主频率信息吗？请使用max函数记录下主频信息，

可以使用help来查询freqz与max来获得更多帮助信息。

2.因为干扰为正弦波，缺穗它应该具有如下形式：n(t) = Acos(wt +f) 其中w =800p 振幅A 相位f未知。取样后（这里不理解，我认为是一种离散化的数值表示），

离散时间的干扰信号是n[k ] = Acos(wkTs +f) 式中：Ts =1/ 22000 经过三角恒等变换 n[k ] =Acos(f)cos(wkTs ) -A sin(f)sin(wkTs )

=an'1 [k ] +bn'2 [k ]其中n'1 [k ] = cos(wkTs ) and n'2 [k ] = sin(wkTs ) 可以被当做求扮渗解的原始已知信息，a =Acos(f) and b=-A sin(f)

是应该被考虑的未知权重项，利用最小二乘法对a，b进行估计：（这里是不少打字了）

（d）得到a，b后，编写matlab代码恢复原始文件依据sˆ[k ] =r[k ] -aˆ n'1 [k ] -bˆ n'2 [k ]表达式。

#include<graphics.h>

#include<stdlib.h>档颂祥

#include<stdio.h>

#include<math.h>

main()

{

int graphdrv=VGA

int graphmode=VGAHI

int i

float k=0.75,dt=0.05,x0,y0,x,y,vx,vy,ax,ay,r1,r2

initgraph(&graphdrv,&graphmode,"\\tc\\bgi")

line(400,0,400,480)line(0,240,640,240)

fillellipse(400,160,20,20)settextstyle(1,0,0)

fillellipse(400,320,20,20)settextstyle(1,0,0)

outtextxy(180,0,"SCATTERED")

outtextxy(180,25,"ALPHA RAY")

settextstyle(3,0,0)

for(i=0i<行搏11i++)

{

y0=80+20*i

y=y0

x=x0=2

vx=1.5

vy=0

r1=sqrt(pow(x/10-40,2)+pow(y/樱颂10-16,2))

r2=sqrt(pow(x/10-40,2)+pow(y/10-32,2))

ax=k*((x/10-40)/pow(r1,3)+(x/10-40)/pow(r2,3))

ay=k*((x/10-16)/pow(r1,3)+(y/10-32)/pow(r2,3))

vx=vx+0.5*ax*dt

vy=vy+0.5*ay*dt

loop:x=x+vx*dt

y=y+vy*dt

r1=sqrt(pow(x/10-40,2)+pow(y/10-16,2))

r2=sqrt(pow(x/10-40,2)+pow(y/10-32,2))

ax=k*((x/10-40)/pow(r1,3)+(x/10-40)/pow(r2,3))

ay=k*((y/10-16)/pow(r1,3)+(y/10-32)/pow(r2,3))

vx=vx+0.5*ax*dt

vy=vy+0.5*ay*dt

if(x>639){goto end}

if(x<2){goto end}

if(y<2){goto end}

line(x0,y0,x,y)

line(x0,480-y0,x,480-y)

x0=x

y0=y

goto loop

end:continue

}

while(!kbhit())

closegraph()

}

---