[ViewVC] Diff of: JSOC/proj/globalhs/sosh/soshdata.py

Comparing proj/globalhs/sosh/soshdata.py (file contents):
Revision 1.3 by tplarson, Mon Mar 4 17:33:22 2019 UTC vs.
Revision 1.4 by tplarson, Mon Apr 8 16:52:19 2019 UTC

+import numpy as np
+import os
-+
+import re
+import drms
+from astropy.io import fits
+from urllib.request import urlretrieve
+import soundfile as sf
-+
+from scipy import fftpack
+import warnings
+warnings.simplefilter(action='ignore', category=FutureWarning)
+      print("Invalid number, try again.")
+      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'))
+  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)):
+  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
-+
-+

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing proj/globalhs/sosh/soshdata.py (file contents): Revision 1.3 by tplarson, Mon Mar 4 17:33:22 2019 UTC vs. Revision 1.4 by tplarson, Mon Apr 8 16:52:19 2019 UTC

Diff Legend

Comparing proj/globalhs/sosh/soshdata.py (file contents):
Revision 1.3 by tplarson, Mon Mar 4 17:33:22 2019 UTC vs.
Revision 1.4 by tplarson, Mon Apr 8 16:52:19 2019 UTC