IT博客汇
  • 首页
  • 精华
  • 技术
  • 设计
  • 资讯
  • 扯淡
  • 权利声明
  • 登录 注册

    当教室突然安静·二——衰减与“静默波”

    BaoPinsUi啊发表于 2024-02-19 17:22:32
    love 0

    前情提要:

    https://zhuanlan.zhihu.com/p/679862117

    在上篇文章中,我们完成了模型的建立并使用python进行了模拟。然而,受制于python的动态性导致的速度低下,我们只能模拟 S(f_1,f_2…,f_N)=\frac{f_1,f_2,…,f_N}{N} 的简单情况,而不能模拟真正的传播过程中声音的衰减。于是我使用C++重写了模型。在当前的模型中,传播函数S_{ij}=归一系数*\mathop{\sum}\limits_{m,n}{\frac{f_{mn}}{(0.1r+1)^3}},其中r=\sqrt{(m-i)^2+(n-j)^2},若i=m或j=n则该项为零以避免自身影响。其中i与j是接收者的坐标,m与n是发送者的坐标。 主要代码如下:

    #include <bits/stdc++.h>
    using namespace std;
    int length=10;//人阵规模,不得超出10000*10000
    int width=10;
    //下为f二三项
    double A(double a,double b)
    {
    	double y=log1p(a)/log(2)*0.5*(1+tanh(5*(b+1.5)));
    	return y;
    }
    //下为随机函数
    double g(double t)
    {
    	double y=1+tanh(7*(sin((sqrt(3))*0.7*t)+sin((sqrt(5))*0.7*t)+0.332*sin((sqrt(16))*0.7*t)+sin((sqrt(14))*0.7*t)+1.02*sin((sqrt(2.5803)*0.7*t))));
    	return y;
    }
    //这个g2是用来错开各个f的随机初始值,为了省事直接把g里面的部分拿出来了
    double g2(double t)
    {
    	double y=(sin((sqrt(3))*0.7*t)+sin((sqrt(5))*0.7*t)+0.332*sin((sqrt(16))*0.7*t)+sin((sqrt(14))*0.7*t)+1.02*sin((sqrt(2.5803)*0.7*t)));
    	return y;
    }
    //递推步进函数
    double Nextf(int step,vector<double> ramvec,double deltat)
    {
    	double t=0.05*step;
    	double fl=(ramvec[(step+19)%20]+ramvec[(step+18)%20]+ramvec[(step+17)%20]+ramvec[(step+16)%20]+ramvec[(step+15)%20]+ramvec[(step+14)%20])/6;
    	double fv=(((ramvec[(step+19)%20]+ramvec[(step+18)%20]+ramvec[(step+17)%20]+ramvec[(step+16)%20]+ramvec[(step+15)%20])/5)-((ramvec[(step+9)%20]+ramvec[(step+8)%20]+ramvec[(step+7)%20]+ramvec[(step+6)%20]+ramvec[(step+5)%20])/5))*2;
    	double ft=g(t-deltat)*A(fl,fv)+0.001*g(t-deltat);
    	return ft;
    }
    
    double fdata[10000][10000];//传出响度定义
    //卷积核(迫真)就是模型里的S(……),决定声音如何传播,衰减
    double core(int i,int j)//接收坐标
    {
    	double fsum=0;
    	double f1sum=0;
    	double daor=0;
    	for(int m=0;m<width;m++)//遍历发射坐标
    	{
    		for(int n=0;n<length;n++)
    		{
    			double r=sqrt(pow((m-i),2)+pow((n-j),2));//两人距离(平直空间,不整花活)
    			if(r==0)//排除自身影响
    			{
    				daor=0;	
    			}
    			else
    			{
    				daor=(pow((0.1*r+1),-3));
    			}
    			fsum=fsum+(fdata[m][n])*daor;//求和
    			f1sum=f1sum+daor;//归一化倒系数
    		}
    	}
    	double loud=fsum/f1sum;
    	return loud;
    }
    //——————————————————————————施法前摇结束——————————————————————————————————————————————■■■■■■■■■■■■■■
    int main()
    {
    	//输出
    	freopen("outputtest.txt","w",stdout);
    	//传出响度初始化
    	for(int i=0;i<width;i++)
    	{
    		for(int j=0;j<length;j++)
    		{
    			fdata[i][j]=1;
    		}
    	}
    	
    	vector<double> fsave[width][length];//历史接受响度初始化
    	vector<double> example(20,1);
    	for(int i=0;i<width;i++)
    	{
    		for(int j=0;j<=length;j++)
    		{
    			fsave[i][j]=example;
    		}
    	}
    	
    
    	//数据记录用量(暂时没用到
    //	vector<double> tdata={};//突变时间记录
    //	int count=0;//突变次数
    //	int test=0;
    //	int last=0;
    	long long int step;//这里step指第几步而非步长,步长固定在0.05
    	//———————————————————————————————————递推循环————————————————醒目的分界线————————————————————————————■■■■■■■■■■———— 
    	for(step=0;step<=(20*1000);step=step+1)//迭代循环,1000s用15s
    	{
    		double t=0.05*step;
    		//这里放迭代主体
    		//这里放卷积
    		for(int i=0;i<width;i++)//开卷!
    		{
    			for(int j=0;j<length;j++)
    			{
    				fsave[i][j][step%20]=core(i,j);
    			}
    		}
    		for(int i=0;i<width;i++)//反馈计算
    		{
    			for(int j=0;j<length;j++)
    				{
    					fdata[i][j]=Nextf(step,fsave[i][j],(114*g2(i)+1919*g2(j)));//最后一个自变量用来错开随机
    				}
    		}
    
    		//突变判断(没用上)
    //		if(fall>=0.1)
    //		{
    //			test=0;
    //		}
    //		else
    //		{
    //			test=1;
    //		}
    //		
    //		if(test-last==1)
    //		{
    //			count=count+1;
    //			tdata.push_back(t);
    //		}
    //		last=test;
    		
    		for (int i = 0; i < width; ++i) {//把每个step的矩阵拍扁后输出,这里代码排版莫名其妙炸了
        	    for (int j = 0; j < length; ++j) {
       	         cout<<fdata[i][j]<<" ";
      	      }
      	  }
      	  cout<<""<<endl;
    	}//循环结束——————————————————————————————————————————————■———————很醒目的分界线—————————■■■■■■■■■■■■■■■■■■■■■	
    
    	
    return 0;	
    }

    随后我们把输出的数据丢给py进行可视化[[[[[代码在文末!!!!!!!!!!!]]]]]后便得到了结果。接下来我们先运行一个简单的模型,10*10的人阵,结果如下:

    参数n=1.0,开放边界条件。以防忘记,n是描述警觉阈值的参数,n越小人越警觉。

    结果与上篇文章中的声音曲线相一致:在突变发生前房间的左上角出现了一个大小在3*4左右的高音量区,随后一同静默,并引发了一次突变。

    左图,“偶然波峰”———右图,“偶然波谷”

    然而,这次我们观察到了更有趣的现象。与无音量衰减的模型不同,考虑衰减的模型中静默的产生是有一个传递的过程的,我称其为“静默波”。

    n=1时,波速约为每秒50倍人间距

    为了更显著的观察这一现象,我又进行了参数n=0.5和n=0.35(由于模拟的时间复杂度为O(N²)(这里N是人数),我不得不过度的调小n来促使突变现象在尽可能短的时间内发生)下规模100*100——一万人的模型的模拟,结果如下

    n=0.5,注意右下和右上相继出现的两次静默,右上的形态接近上篇文章中的“不完全静默”
    n=0.35,右侧和左侧两列波相遇后向垂直方向传播,而没有发生相互穿过的现象,这一点与神经冲动很相似,波前后方的人处于一小段静默期中,不支持波穿过后继续传递。图中波速约每秒125倍人距

    现在我们有了更多有趣的现象。在n=0.5的模型中,人阵的右下,右上角分别出现了两次局域的静默,然而都没有传播开来。仔细观察可以发现右下第一个静默波的传播范围和起始波谷剧烈程度都大于第二个静默波;n=0.35的模型中静默波得以持续的传播出去,并且在相遇后展现出来与神经冲动的传导相似的“相遇抵消”性,并且波速几乎三倍与n=1.0的情况。两次模拟中波源都位于边框和角落,这很好理解,位于角落的人身边人更少,人少则波谷相遇概率更高。

    n=0.35,静默波波前用红线标注

    由此我们可以推断:1.静默波的出现频率正相关于敏感程度,正相关与声音衰减速度;2.静默波的传播范围正相关于“振幅”,正相关与人的警觉程度,负相关于声音衰减程度;3.静默波的波速正相关于警觉程度,负相关与声音本身的衰减程度,当声音不衰减时静默波传播速度为无穷大。这是很好定性的解释的:1.人越警觉,越容易静默,而声音快速衰减使得人主要感知的声源越少,波谷就越容易叠加;2.人对静默的差距有0.3秒左右的延迟,当n偏大时前一波人静默与后一波人静默的时差会抹平突变的剧烈程度,限制静默波的范围,而提高波源强调会延缓其发生,表现为传播距离增加;3.人越警觉,声音越不衰减,人就越能敏感的感觉到远处声音的变化,每一步模拟静默波传播的范围就更广,波速就更快。可惜受限于算力,我无法定量的测量它们之间互相影响的具体形式。

    最后,请欣赏: 准 一 维 人 长 链(5*200的人方阵,原则上在这种狭小空间声音不该衰减,就当是露天会议罢((

    https://www.zhihu.com/video/1738355168486293504

    (gif太大放不上来,为了比例正常我把5的那一边拉宽了亿点,我期望看到神经冲动一样波对消的场景,可惜没遇到)。

    最后放上python可视化部分的代码(py可视化费时比c艹模拟数据慢多了)

    import os
    import math
    import numpy as np
    import matplotlib.pyplot as plt
    current_dir = os.getcwd()
    
    folder_path = "/storage/emulated/0/Documents/Pydroid3"
    file_name = "long.txt"
     
    full_path = os.path.join(current_dir, folder_path, file_name)
    
    with open(full_path, 'r') as file:
    	content = file.read().splitlines() 
    	list = [[float(x) for x in line.split()]for line in content] 
    fig=plt.figure(figsize=(10,10))
    fig,ax = plt.subplots()
    for step in range(0,1200):#把被拍扁的矩阵立起来
    	for i in range(0,1000):
    		a=math.tanh(list[step][i])
    		rect=plt.Rectangle((i%200,(i//200)),1,1,alpha=a)
    		ax.add_patch(rect)
    		plt.axis([0, 200, 0, 5])
    	
    	plt.savefig('/storage/emulated/0/Documents/Pydroid3/longaa/'+str(step)+'.png')
    	ax.cla()
    	
    
    	

    好了,这就是当教室突然安静系列的全部主要内容,感谢各位观看。以后若有点子,我会随缘更新。



    来源:知乎 www.zhihu.com
    作者:BaoPinsUi啊

    【知乎日报】千万用户的选择,做朋友圈里的新鲜事分享大牛。 点击下载


沪ICP备19023445号-2号
友情链接