2017年12月11日月曜日

人間の耳で音の位相は聞き取れるのか?

「人間の耳で音の位相は聞き取れるのか?」問題のテスト。音痴な自分でも結構分かる気はする...(特にphase=0と0.5の差は相当顕著)。
import ddf.minim.*;    
import ddf.minim.ugens.*;    

Minim minim;
AudioOutput out;    
Oscil sin1, sin2;
float f=440;
float A=1.0;
float trig=0.5;
void setup() {

  size(512, 200);
  smooth();
  minim = new Minim(this);
  out = minim.getLineOut(Minim.MONO, 1024);

  sin1=new Oscil(f, A, Waves.SINE);
  sin2=new Oscil(f*2, A*0.5, Waves.SINE);
  sin1.setPhase(0);
  sin2.setPhase(0);

  sin1.patch(out);
  sin2.patch(out);
}


void draw() {
  background(255);
  stroke(0);
  float x=0, Amp_old=999;
  int flag=0;
  for (int i = 0; i < out.bufferSize()-1; i++) {
    float Amp=out.left.get(i);
    if ((Amp_old=trig)||flag==1) {
      point(x, height/2 + Amp*height/4);
      flag=1;
      x+=2;
    }
    Amp_old=Amp;
  }
}

void mouseMoved() {
  float phase=map(mouseX, 0, width, 0, 1);
  sin2.setPhase(phase);
}


void stop() {
  out.close();
  minim.stop();
  super.stop();
}

2017年12月7日木曜日

Pyorbitalで方位仰角、Rangeの計算

Pyorbitalで方位仰角、Rangeの計算。こんなのがさっくりできるのは、とてもありがたい...。どのくらい信用できるのかはわからないけど。
from pyorbital import orbital #これを追加すると動く(外すと動かない)のは何故
from pyorbital import orbital
import pyorbital
import numpy
import datetime

#Sat info
satname='ITUPSAT 1'
f=437.325
tlef='TLE.txt'

#Position Info
#静岡大学
lon=138.431704
lat=34.965720
alt=90

#Obs Info
start_jst=datetime.datetime(2018,02,16,15,25,00)
end_jst=datetime.datetime(2018,02,16,15,45,00)

#Conf
c=299792.458 # km/s
dt=datetime.timedelta(0,30)

def jst2utc(time) :
    return time-datetime.timedelta(0,9.*60*60,0)

def obsf_ICR6(obsf) :
    f0=obsf*1000.
    f1=int(obsf*100.)*10
    df=f0-f1
    if df <2.5 :
        return f1/1000.
    elif df <7.5 : 
        return (f1+5.)/1000. 
    else :
        return (f1+10.)/1000. 
    
tle=pyorbital.tlefile.read(satname,tlef)
orb=pyorbital.orbital.Orbital(satname,tlef)

print "TLE Epoch:",",",str(datetime.datetime(2000+int(tle.epoch_year),1,1)+datetime.timedelta(days=tle.epoch_day)).split(' ')[0]
print "Date:",",",str(start_jst).split(' ')[0]
print "Lon:",",",lon
print "Lat:",",",lat
print "Alt:",",",alt
print ""

jst=start_jst
while jst < end_jst :

    now=jst2utc(jst)
    obs_pos=numpy.array(pyorbital.astronomy.observer_position(now,lon,lat,alt)[0])
    sat_pos,sat_vel=orb.get_position(now,False)

    range_dis=numpy.linalg.norm(sat_pos-obs_pos)
    range_vec=(sat_pos-obs_pos)/range_dis
    rang_rate=numpy.dot(sat_vel,range_vec)

    Az,El=orb.get_observer_look(now,lon,lat,alt)
    obsf=f*(1-rang_rate/c)
    
    print str(jst).split(' ')[1],",",Az,",",El,",",obsf_ICR6(obsf),",",obsf
    jst+=dt

2017年6月26日月曜日

大学入試問題をモンテカルロで解く

情報処理の課題として、「大学入試問題をモンテカルロで解く」というのは面白いかと思い、こちらの問題について作成。しかし、これなら素直に数え上げた方が遥かに速い(汗)。


int a,b,c;
float n1=0;
float n2=0;
void setup(){
  size(1000,20);
}
void draw(){
   a=int(random(1,7));
   b=int(random(1,7));
   c=int(random(1,7));
   
   if(log(a+b)/log(0.25)>log(c)/log(0.5)) n1+=1;
   if((pow(2,a)+pow(2,b)+pow(2,c))%3==0) n2+=1;
   background(255);
   fill(0);
   text(str(frameCount),width*0.1,height/2);
   text(str(n1/frameCount),width*0.5,height/2);
   text(str(n2/frameCount),width*0.8,height/2);
}

数えた。コンピュータ便利。
int a, b, c;
float n1=0;
float n2=0;
float n=0;
size(1000, 20);
for (a=1; a<7; a+=1) {
  for (b=1; b<7; b+=1) {
    for (c=1; c<7; c+=1) {
      if (log(a+b)/log(0.25)>log(c)/log(0.5)) n1+=1;
      if ((pow(2, a)+pow(2, b)+pow(2, c))%3==0) n2+=1;
      n+=1;
    }
  }
}
background(255);
fill(0);
text(str(n), width*0.1, height/2);
text(str(n1/n), width*0.5, height/2);
text(str(n2/n), width*0.8, height/2);

2017年5月2日火曜日

画像データの値を取得

こちらを参考に作成。面白い
PImage myPhoto;
int myPhotoWidth, myPhotoHeight;

void setup() {
  size(500,750);
  //Thanks for http://gahag.net/004881-high-school-girl/
  myPhoto = loadImage("Test.jpg"); 
  myPhotoWidth = myPhoto.width;
  myPhotoHeight = myPhoto.height;
  noStroke();
  smooth();
  background(0);
  myPhoto.loadPixels();
}

void draw() {
  background(0,0,100);

  int x;
  int y; 
  for(x=0;x<myPhoto.width;x+=5){
  for(y=0;y<myPhoto.height;y+=5){
  
  strokeWeight(random(1,3));
  float f=random(0.5,1);
  stroke(red(myPhoto.pixels[y*width + x])*f,green(myPhoto.pixels[y*width + x])*f,blue(myPhoto.pixels[y*width + x])*f);

  float r0=random(5, 15);
  float r1=random(5, 15);
  float th0=random(0, 2*PI);
  float th1=random(0, 2*PI);
  float dx=random(-5, 5);
  float dy=random(-5, 5);

  line(r0*cos(th0)+x+dx, r0*sin(th0)+y+dy, r1*cos(th1)+x+dx, r1*sin(th1)+y+dy);
   }}
}

2017年4月30日日曜日

色々統合した練習

もう少し賢く書けそう。こちら
PImage img;
int R=200, G=200, B=200;

void setup() {
  size( 800, 800 );
  background( R, G, B );
  //Thanks for Yume Yume Iro TOWN http://www.poipoi.com/yakko/cgi-bin/sb/log/eid1737.html
  img = loadImage("siruetto20080205a.png");
  stroke(255, 0.0);
  noCursor();
}
void draw() {
  background(R, G, B);
  float x=mouseX-280;
  float y=mouseY-500;


  for (int i=0; i<10; i++) {
    fill(R, G, B, 30);
    noStroke();
    rect(0, 0, width, height);
    float xi=map(i, 0, 10, width/2, x);
    float yi=map(i, 0, 10, height/2, y);
    float f=map(i, 0, 10, 0.1, 1);
    image(img, xi, yi, 560*f, 1000*f);
  }
}
void mouseMoved() {
  int a=int(random(1, 4));
  int pm;
  if (R+G+B<50) {
    pm=1;
  } else if (random(0, 1.0)<0.5) {
    pm=-1;
  } else {
    pm=1;
  }

  if (a==1) {
    if (R==200) {
      R=199;
    } else if (R==1) {
      R=2;
    } else {
      R=R+pm;
    }
  }
  if (a==2) {
    if (G==200) {
      G=199;
    } else if (G==1) {
      G=2;
    } else {
      G=G+pm;
    }
  } 
  if (a==3) {
     if (B==200) {
      B=199;
    } else if (B==1) {
      B=2;
    } else {
      B=B+pm;
    }
  } 
}

2017年4月29日土曜日

マウスの使い方練習

マウスの使い方練習。IRHH 2
PImage img;

void setup() {
  size( 600, 600 );
  background( 0, 0, 255 );
  //Thanks for Yume Yume Iro TOWN http://www.poipoi.com/yakko/cgi-bin/sb/log/eid1737.html
  img = loadImage("siruetto20080205a.png");
  stroke(255, 0.0);
}
void draw() {
  background(0, 0, 255);
  image(img, mouseX-140, mouseY-250, 280, 500);
  //image(img, mouseX-120, mouseY-80, 280, 500);

  beginShape();
  for (int i=0; i<300; i++) {
    stroke(random(0, 255), 0, 0);
    strokeWeight( random(1, 3) );

    float r0=random(0, 50);
    float r1=random(0, 50);
    float th0=random(0, 2*PI);
    float th1=random(0, 2*PI);
    line(r0*cos(th0)+mouseX-25, r0*sin(th0)+mouseY-170, r1*cos(th1)+mouseX-25, r1*sin(th1)+mouseY-170);
  }
}

ペジエの使い方の練習

ペジエの使い方の練習。こちら
PImage img;


void setup() {
  size( 800, 800 );
  background( 255 );
  //Thanks for Yume Yume Iro TOWN http://www.poipoi.com/yakko/cgi-bin/sb/log/eid1737.html
  img = loadImage("siruetto20080205a.png");
  strokeWeight( 6 );
}

void draw() {
  if (frameCount%3==0) {
    float f=random(0.1, 0.8);
    image(img, random(0, 500), random(0, 500), 560*f, 1000*f);
  }
  stroke(  random(255), random(255), random(255) );
  fill( random(255), random(255), random(255));


  beginShape();
  int x0=radomxy();
  int y0=radomxy();

  vertex( x0, y0 );
  bezierVertex( radomxy(), radomxy(), radomxy(), radomxy(), radomxy(), radomxy() );
  int i=0;
  while (i==0) {
    if (random(6)>3) {
      i=1;
    }
    bezierVertex( radomxy(), radomxy(), radomxy(), radomxy(), radomxy(), radomxy() );
  }
  bezierVertex( radomxy(), radomxy(), radomxy(), radomxy(), x0, y0 );
  endShape();
}

void mousePressed() {
  //background( 255 );
  delay(1000);
}

int radomxy() {
  return int(random(-100, height+100));
}

2017年4月26日水曜日

なにかとむなしい

こちら

int x, y, vx, vy;

void setup() {
  smooth();
  size(800, 800);
  x=width/2;
  y=height/2;
  background(255);
  int v=20;
  vx=int(random(-1*v, v));
  int a;
  if (random(1)>0.5) {
    a=1;
  } else {
    a=-1;
  }

  vy=a*int(sqrt(pow(v, 2)-pow(vx, 2)));
}

void draw() {
  float a=random(1, 4);
  int xnew=x+int(a*vx);
  int ynew=y+int(a*vy);
  if (xnew<0) {
    xnew=0;
    vx=-1*vx;
    vy=int(vy*random(0.9, 1.1));
  }
  if (xnew>width) {
    xnew=width;
    vx=-1*vx;
    vy=int(vy*random(0.9, 1.1));
  }
  if (ynew<0) {
    ynew=0;
    vy=-1*vy;
    vx=int(vx*random(0.9, 1.1));
  }
  if (ynew>height) {
    ynew=height;
    vy=-1*vy;
    vx=int(vx*random(0.8, 1.2));
  }
  if (vx*vx<30) {
    vx=int(random(-20, 20));
  }
  if (vy*vy<30) {
    vy=int(random(-20, 20));
  }  
  stroke(random(255), random(255), random(255));
  strokeWeight(random(30));
  line(x, y, xnew, ynew);
  x=xnew;
  y=ynew;
  if (frameCount%10==0) {
    fill(255, 20);
    noStroke();
    rect(0, 0, width, height);
  }
}

2017年4月19日水曜日

月の自転・公転周期デモ

月の自転・公転周期の説明デモ。動いているのはこちら。意図せずもドラoもんカラー。
float theta=0;
int n=0;

void setup()
{
  size(1000, 1000);
  background(255);
  noLoop();
}

void draw()
{
  float r=width/3;
  float re=width/8;
  float rm=width/10;

  background(255);
  stroke(100);
  fill(0, 0, 255);
  ellipse(width/2, height/2, re, re);
  stroke(0);

  fill(255, 255, 0);

  ellipse(width/2+r*cos(theta), height/2+r*sin(theta), rm, rm);
  fill(255, 0, 0);

  ellipse(width/2+r*cos(theta)+0.35*rm*cos(PI+1*theta), height/2+r*sin(theta)+0.35*rm*sin(PI+1*theta), rm*0.2, rm*0.2);
  theta+=PI/200;
}

void mousePressed() {
  if (n==0) {
    loop();
    n=1;
  } else {
    noLoop();
    n=0;
  }
}  

2017年4月14日金曜日

マル9円周率その3

意味無しモンテカルロ。こちら
float n=1;
float s=0;

void setup() {
  size(300, 350);
  background(255);  
  noFill();
  stroke(0, 2555, 0);
  ellipse(0, 0, 600, 600);
}

void draw() {

  float x=random(0.0, 1.0);
  float y=random(0.0, 1.0);
  if (pow(x, 2)+pow(y, 2)<1.0) {
    stroke(255, 0, 0);
    s+=1;
  } else {
    stroke(0, 0, 255);
  }
  ellipse(x*300, y*300, 1, 1);
  stroke(255);
  fill(255);
  rect(0, 310, 300, 50);
  fill(0);
  text(str(int(n)), 5, 330);  
  text(str(s/n*4.0), 150, 330);
  n++;
}

void mousePressed() {
  n=1;
  s=0;
  background(255);
  noFill();
  stroke(0, 2555, 0);
  ellipse(0, 0, 600, 600);
}  

マル9な円周率その2

その2。動いているのはこちら
float n=1;

void setup(){
  size(300,350);
  background(255);  
}

void draw(){
  float s=0;
  background(255);
  fill(255);
  noFill();
  stroke(255,0,0);
  ellipse(0,0,600,600);
  float x0=0;
  float y0=1.0;


  for(int m=1;m<n+1;m++){
    float x=m/n;
    float y=sqrt(1-pow((x),2.0));
    stroke(0);
    line(x*300,y*300,x0*300,y0*300);
    s=s+sqrt(pow((x-x0),2.0)+pow((y-y0),2.0));
    x0=x;
    y0=y;
}
  n+=int(n/10+1);
 if(n<300){
   delay(int(1000/sqrt(n)));
  }
  fill(0);
  textSize(15);
  text(str(int(n)),5,330);  
  text(str(s*2.0),150,330);
  if(n>10000000){
    noLoop();
  }    
}

void mousePressed(){
  n=1;
  loop();
}  

2017年4月13日木曜日

マル9な円周率計算デモ

頭の悪い円周率計算デモ。Processing.jsで動かすとこちら
float n=1;

void setup(){
  size(300,350);
  background(255);  
}

void draw(){
  float s=0;
  background(255);
  fill(255);
  for(int m=0;m<n;m++){
    float x=m/n;
    float f=sqrt(1-pow((x),2.0));
    stroke(0);
    rect(x*300,0,1/n*300,f*300);
    noFill();
    stroke(255,0,0);
    ellipse(0,0,600,600);

    s=s+f*1/n;
  }
  n+=int(n/10+1);
 if(n<300){
    delay(int(1000/sqrt(n)));
  }
  fill(0);
  textSize(15);
  text(str(int(n)),5,330);  
  text(str(s*4.0),150,330);
  if(n>10000000){
    noLoop();
  }    
}

void mousePressed(){
  n=1;
  loop();
}  

2017年4月7日金曜日

地球の自転に対するアリストテレス or ニュートン的描像

天文学の授業で「地動説に対する反論」を説明する為にProcessingで作成したアニメーション。こういう教材を簡単に作れるので、本当はどこかの授業で学生さんにProcessingを触ってもらいたいところ。Processing.jsで動かしたアリストテレス or ニュートン的な描像はそれぞれ、こちらこちら
float glevel=0.8;

float x;

float xb, yb, vy=0, vxw;

float headsize;

int t=0;

int ball_start=150;

float g;

float[] xbo=new float[500];

float[] ybo=new float[500];



void setup() {

  x=100;

  size(1100, 700);

  colorMode(RGB, 256);

  background(255, 255, 255);

  headsize=width*0.03;

  frameRate(50);

  xb=0;

  yb=height*glevel*0.15;

  g=0.05;

  vxw=width*0.002;



  for (int i=0; i<500 i="" p="">    xbo[i]=0;

    ybo[i]=0;

  }

}

void draw() {

  background(255, 255, 255);

  //smooth() ;

  //World

  stroke(0, 0, 0);

  strokeWeight(1);

  line(0, height*glevel, width, height*glevel);

  noStroke();

  //point

  fill(255, 0, 0);

  ellipse(x, height*glevel, width*0.01, width*0.01);

  //man

  fill(0, 0, 0);



  //head

  ellipse(x+width*0.1, height*glevel*0.71, headsize, headsize);

  //body

  rect(x+width*(0.1-0.025), height*glevel*0.75, width*(0.05), height*0.1);

  //leg

  rect(x+width*(0.1-0.025), height*glevel*0.89, width*(0.02), height*0.085);

  rect(x+width*(0.1+0.005), height*glevel*0.89, width*(0.02), height*0.085);

  //arm

  rect(x+width*(0.1-0.045), height*glevel*0.75, width*(0.015), height*0.08);

  rect(x+width*(0.1+0.03), height*glevel*0.75, width*(0.015), height*0.08);



  //Bar

  strokeWeight(10);

  stroke(183, 65, 14);

  line(x+width*0.2, height*glevel, x+width*0.2, height*glevel*0.1);

  line(x, height*glevel*0.1, x+width*0.2, height*glevel*0.1);



  //Ball Shadow

  for (int i=1; i    if ( i%20==0) {

      stroke(150);

      strokeWeight(1);

      noFill();

      ellipse(xbo[i], ybo[i], width*0.03, width*0.03);

    }

  }



  //Ball

  xbo[t]=xb;

  ybo[t]=yb;

  if (t>ball_start & yb<=height*glevel) {

    vy=vy+g;

    yb=yb+vy;

    //下の1文を入れるか入れないかが両者の差。

    xb=x;

  } else {

    xb=x;

  }



  if (yb>=height*glevel) {

    vy=0;

    yb=height*glevel;

    xb=xb+vxw;

  }





  noStroke();

  fill(0, 0, 255);

  ellipse(xb, yb, width*0.03, width*0.03);

  x=x+vxw;

  print(t, "\n");

  if (t<499 p="">    t++;

  }

}

2017年4月5日水曜日

Processingで矩形波のフーリエ変換の可視・音声化デモ

こちらを参考にさせていただきました。学生さんの卒論に向けたProcessingでのフーリエ合成(というかMinimの使い方)練習。Processing、本当に楽しい...。

import ddf.minim.*;   
import ddf.minim.signals.*;   


Minim minim;
AudioOutput out;  
SineWave[] sine;
float waveH = 50;  
int n=30;



void setup()
{

  size(512, 200);
  smooth();
  minim = new Minim(this);
  sine=new SineWave[n];

  out = minim.getLineOut(Minim.STEREO);

  int i;

  for(i=0;i < n;i++){
    int m=i+1;
    sine[i] = new SineWave(440*m, 2*(1-cos(m*PI))/(m*PI*440)*100, out.sampleRate());
  }



  for(i=0;i < n;i++){
  out.addSignal(sine[i]);
  }

}



void draw()
{
  background(255);
  stroke(0);
  for(int i = 0; i < out.bufferSize()-1; i++)

  {

    point(i, 50 + out.left.get(i)*waveH); 
    point(i, 150 + out.right.get(i)*waveH);

  }

}

void mouseMoved() {
  int t = int(map(mouseX, 0, width, 1, n));

  for(int i=0;i < t;i++){
    int m=i+1;
    sine[i].setAmp(2*(1-cos(m*PI))/(m*PI*440.)*100); 

  }
  for(int i=t;i < n;i++){
    sine[i].setAmp(0); 
  }
}



void stop()
{
  out.close();
  minim.stop();
  super.stop();
}