(file) Return to soshdata.py CVS log (file) (dir) Up to [Development] / JSOC / proj / globalhs / sosh

Diff for /JSOC/proj/globalhs/sosh/soshdata.py between version 1.3 and 1.4

version 1.3, 2019/03/04 17:33:22 version 1.4, 2019/04/08 16:52:19
Line 1 
Line 1 
 import numpy as np import numpy as np
 import os import os
   import re
 import drms import drms
 from astropy.io import fits from astropy.io import fits
 from urllib.request import urlretrieve from urllib.request import urlretrieve
 import soundfile as sf import soundfile as sf
   from scipy import fftpack
  
 import warnings import warnings
 warnings.simplefilter(action='ignore', category=FutureWarning) warnings.simplefilter(action='ignore', category=FutureWarning)
Line 168  def getwavfiles(instrument=None, day=Non
Line 170  def getwavfiles(instrument=None, day=Non
       print("Invalid number, try again.")       print("Invalid number, try again.")
       return None       return None
  
   wblank = '%s.%id_l=%i_m=%i_data%s.wav'    wblank = '%s.%id_72d.l=%i_m=%i_data%s.wav'
   wfile=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'r'))   wfile=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'r'))
   wfile2=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'i'))   wfile2=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'i'))
   if os.path.exists(wfile) and (m==0 or os.path.exists(wfile2)):   if os.path.exists(wfile) and (m==0 or os.path.exists(wfile2)):
Line 293  def apols(lin, amaxin):
Line 295  def apols(lin, amaxin):
   return pols   return pols
  
  
   def splitwavfile(wavfile=None, splitfactor=2, splitboth=True):
   
     wblank = '%s.%id_%id.l=%i_m=%i_data%s.wav'
     outlist=[]
   
     if (wavfile == None):
       inst=input("Enter name of instrument: ")
       inst=inst.lower()
       if (inst != 'mdi' and inst != 'hmi'):
         print("Invalid instrument, try again.")
         return None
       daystr=input("Enter day number: ")
       try:
         firstday=int(daystr)
       except:
         print("Invalid number, try again.")
         return None
       ndstr=input("Enter number of days: ")
       try:
         ndaysin=int(ndstr)
       except:
         print("Invalid number, try again.")
         return None
       lstr=input("Enter spherical harmonic degree: ")
       try:
         l=int(lstr)
         if (l < 0):
           print("Degree l must be nonnegative, try again.")
           return None
       except:
         print("Invalid number, try again.")
         return None
       mstr=input("Enter azimuthal order: ")
       try:
         m=int(mstr)
         m=np.abs(m)
         if (m > l):
           print("Abs(m) must be less than l, try again.")
           return None
       except:
         print("Invalid number, try again.")
         return None
       ri='r'
       wavfile=os.path.join(datadir,wblank % (inst,firstday,ndaysin,l,m,ri))
     else:
       regmatch = re.search(r"(\S+)\.(\d+)d_(\d+)d\.l=(\d+)_m=(\d+)_data(\S+).wav", wavfile)
       if (regmatch == None):
         print("Invalid file name, try again.")
         return None
       else:
         inst=regmatch.group(1)
         firstday=int(regmatch.group(2))
         ndaysin=int(regmatch.group(3))
         l=int(regmatch.group(4))
         m=int(regmatch.group(5))
         ri=regmatch.group(6)
   
     test= ndaysin % splitfactor
     if (test != 0):
       print("Timeseries not divisible by that factor.")
       return None
   
   #  wfile = os.path.join(datadir,wavfile)
     try:
       x, sr = sf.read(wavfile)
     except:
       print("Invalid file, try again.")
       return None
   
     if (inst == 'mdi'):
       nperday = 1440
     elif (inst == 'hmi'):
       nperday = 1920
   
     ndaysout = ndaysin / splitfactor
     ndt=int(ndaysout*nperday)
     day = firstday
     index=0
     while (day < firstday+ndaysin):
       xout=x[index:index+ndt]
       wfileout = os.path.join(datadir,wblank % (inst,day,ndaysout,l,m,ri))
       sf.write(wfileout, xout, sr, subtype='FLOAT')
       day += ndaysout
       index += ndt
       outlist.append(wfileout)
   
     if (splitboth and m > 0):
       if (ri == 'i'):
         ir = 'r'
       elif (ri == 'r'):
         ir = 'i'
       wfile=wblank % (inst,firstday,ndaysin,l,m,ir)
       wavfile = os.path.join(datadir,wfile)
       try:
         x, sr = sf.read(wavfile)
       except:
         print("Invalid file, try again.")
         return None
       day=firstday
       index=0
       while (day < firstday+ndaysin):
         xout=x[index:index+ndt]
         wfileout = os.path.join(datadir,wblank % (inst,day,ndaysout,l,m,ir))
         sf.write(wfileout, xout, sr, subtype='FLOAT')
         day += ndaysout
         index += ndt
         outlist.append(wfileout)
   
     return outlist
   
   
   def scanspectrum(wavfile=None,downshift=4,firstbin=None,lastbin=None,nbins=200,ndt=10000,binstep=None,ramp=50):
   
     if (binstep == None):
       binstep = int(nbins/2)
   
     wblank = '%s.%id_%id.l=%i_m=%i_data%s.wav'
     if (wavfile == None):
       inst=input("Enter name of instrument: ")
       inst=inst.lower()
       if (inst != 'mdi' and inst != 'hmi'):
         print("Invalid instrument, try again.")
         return None
       daystr=input("Enter day number: ")
       try:
         firstday=int(daystr)
       except:
         print("Invalid number, try again.")
         return None
       ndstr=input("Enter number of days: ")
       try:
         ndaysin=int(ndstr)
       except:
         print("Invalid number, try again.")
         return None
       lstr=input("Enter spherical harmonic degree: ")
       try:
         l=int(lstr)
         if (l < 0):
           print("Degree l must be nonnegative, try again.")
           return None
       except:
         print("Invalid number, try again.")
         return None
       mstr=input("Enter azimuthal order: ")
       try:
         m=int(mstr)
         m=np.abs(m)
         if (m > l):
           print("Abs(m) must be less than l, try again.")
           return None
       except:
         print("Invalid number, try again.")
         return None
       ri='r'
       wavfile=os.path.join(datadir,wblank % (inst,firstday,ndaysin,l,m,ri))
       print("Wav file is %s" % wavfile)
     else:
       regmatch = re.search(r"(\S+)\.(\d+)d_(\d+)d\.l=(\d+)_m=(\d+)_data(\S+).wav", wavfile)
       if (regmatch == None):
         print("Invalid file name, try again.")
         return None
       else:
         inst=regmatch.group(1)
         firstday=int(regmatch.group(2))
         ndaysin=int(regmatch.group(3))
         l=int(regmatch.group(4))
         m=int(regmatch.group(5))
         ri=regmatch.group(6)
   
     try:
       x, sr = sf.read(wavfile)
     except:
       print("Invalid file, try again.")
       return None
   
     z=np.zeros((len(x)),dtype=np.complex_)
     if (ri == 'r'):
       z.real=x
     else:
       z.imag=x
   
     if (m > 0):
       if (ri == 'i'):
         ir = 'r'
       elif (ri == 'r'):
         ir = 'i'
       wfile=wblank % (inst,firstday,ndaysin,l,m,ir)
       wavfile = os.path.join(datadir,wfile)
       try:
         x, sr = sf.read(wavfile)
       except:
         print("Invalid file, try again.")
         return None
     else:
       x=0.0*x
   
     if (ri == 'r'):
       z.imag=x
     else:
       z.real=x
   
     window=np.ones(ndt)
     if (ramp != 0):
       rampsamples=int(sr*ramp/1000)
       if (rampsamples > ndt/2):
         print("Ramp longer than ndt/2, try again.")
         return None
       window[0:rampsamples]=np.linspace(0.0,1.0,rampsamples,endpoint=False)
       window[ndt-1:ndt-rampsamples-1:-1]=np.linspace(0.0,1.0,rampsamples,endpoint=False)
   
     nx=len(x)
     nx2=int(nx/2)
     if (firstbin == None):
       firstbin = int(nx2/6)
     if (lastbin == None):
       lastbin = int(0.9*nx2)
     ftz = fftpack.fft(z)
     freq = fftpack.fftfreq(nx) #*sr
     if (m >= 0):
       mag = np.abs(ftz[0:nx2])
       phase = np.arctan2(ftz[0:nx2].imag,ftz[0:nx2].real)
       f = freq[0:nx2]/downshift
     else:
       mag = np.abs(ftz[nx:nx2:-1])
       phase = np.arctan2(-1*ftz[nx:nx2:-1].imag,ftz[nx:nx2:-1].real)
       f = -1*freq[nx:nx2:-1]/downshift
   
     nstep=int((lastbin-firstbin)/binstep)
     tlen=nstep*ndt
     t=np.linspace(0,ndt-1,ndt)
     xout=np.zeros((tlen))
     bindex=firstbin
     tindex=0
     t0in=0
     while (bindex + nbins <= lastbin):
       for i in range(nbins):
         xout[tindex:tindex+ndt]+=mag[bindex+i]*np.sin(2*np.pi*f[bindex+i]*(t+tindex)+phase[bindex+i])
       xout[tindex:tindex+ndt]*=window
       bindex+=binstep
       tindex+=ndt
   
     mx=np.abs(xout).max()
     xout /= mx
     sf.write('output.wav', xout, sr, subtype='FLOAT')
   
     return xout
   
   
   #def plotspectrum(wavfile):
   #  x, sr = sf.read(wavfile)
   #  ftx = fftpack.fft(x)
   #  freqs = fftpack.fftfreq(len(x)) * sr
   
   

Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
Karen Tian
 Powered by
ViewCVS 0.9.4