Development of a phase velocity
spectral analysis software package(M-transform) for airglow imaging
data and its application on
atmospheric gravity waves studies
Septi Perwitasari1, 4, Takuji Nakamura1, 2, Masaru
Kogure2, 1, Yoshihiro Tomikawa1, 2, Mitsumu K. Ejiri1,2 , Kazuo Shiokawa3
1National Institute of Polar Research, Tokyo, 190-8518, Japan

2Department of Polar Science, SOKENDAI (The Graduate University

for Advanced Studies), Tokyo, 190-8518, Japan

3National Institute of Aeronautics and Space (LAPAN) of Indonesia,Bandung, 40173, Indonesia

4Institute for Space-Earth Environmental Research, NagoyaUniversity, Nagoya, 464-8601, Japan1

Matsuda Transform (M-Transform)
Matsuda et al., 2014
F1 (x, y, t)

F4 (vx, vy)

N

NI’/I
Syowa 11/09/20 23:00

400 kmF3 (k,l,ω)W

3D FFTW

400 km
Log10 PSD [s2.m-2]

2Transform airglow intensity data to power spectrum in horizontal phase velocity domain.

How does it work?
Raw dataGeo. coord

Zero padding

Image preprocessing
Coble et al., 1998
Wavenumber domain

3D FFT

Matsuda et al., 2014
Phase velocity domain

Integration

Wave number tophase speed domain

Integration over frequency

3M-Transform Software

4How to use?
Airglow intensity data at Syowa
station 11/09/20 23:00-06:00
Input=airglow image (400 x 400 km,number of image= 21)
Keywords:
dt=180 (image resolution=3 minutes)
LH= default (5-100 km)
T= default (8-60 min)Vp= default (0-150 m/s)

Airglow data (x,y,t)N

5

Performance Analysis
We examine the performance of this function byrunning several simulations:
qChanging the wave packet size (X) and wave
duration (T) using test data:
1. Case 1: c constant (40 m/s), τ and LH vary (τ=8,15, 30 min).
2. Case 2: c constant, X varies (FWHM=50, 100, 200km).
3. Case 3: c constant (40 m/s), T varies (30, 60, 120min).
q Aliasing effect analysis.
q Zero padding size analysis
6q Noise sensitivity analysis

Simulation Set Up
X=100 km

Ø X, Y, T = Size of wave
packet size in horizontal,meridional and temporal
direction𝒚 = 𝑨 𝒔𝒊𝒏 𝒘𝒂𝒗𝒆𝒏𝒖𝒎𝒙 ∗ 𝑿𝒔 + 𝒘𝒂𝒗𝒆𝒏𝒖𝒎𝒚 ∗ 𝒀𝒔 − 𝝎 ∗ 𝑻𝒔
𝟐𝝅
𝟐𝝅𝐜𝐨𝐬( )𝑳𝑯𝜽
𝟐𝝅
𝟐𝝅𝒘𝒂𝒗𝒆𝒏𝒖𝒎𝒃𝒆𝒓𝒚 =𝐬𝐢𝐧( )𝑳𝑯𝜽𝒘𝒂𝒗𝒆𝒏𝒖𝒎𝒃𝒆𝒓𝒙 =𝝎 = 𝟐𝝅/τ

Ø Xs, Ys, Ts = Size ofairglow image in x (km, y(km) and time (second)Ø Zpx, Zpy, Zpt = Zeropadding size

Case 1: Different wave period and horizontal
wavelength𝒚 = 𝑨 𝒔𝒊𝒏 𝒘𝒂𝒗𝒆𝒏𝒖𝒎𝒙 ∗ 𝑿𝒔 + 𝒘𝒂𝒗𝒆𝒏𝒖𝒎𝒚 ∗ 𝒀𝒔 − 𝝎 ∗ 𝑻𝒔
𝟐𝝅𝟐𝝅𝒘𝒂𝒗𝒆𝒏𝒖𝒎𝒃𝒆𝒓𝒙 =𝐜𝐨𝐬( )𝑳𝑯𝜽
𝟐𝝅𝟐𝝅𝒘𝒂𝒗𝒆𝒏𝒖𝒎𝒃𝒆𝒓𝒚 =𝐬𝐢𝐧( )𝑳𝑯𝜽𝝎 = 𝟐𝝅/τc= 40 m/s
LH= 30 km
τ= 8, 15, 30 min
Ts=120 min
Xs=400 km
dT=1 minΤs_min= 3dT
θ = 225 deg(0=eastward)• The longer the period, the wider thespectrum.
• i.e. The smaller number of cycle
causes the wider spectrum.
• The total power difference is ~30%.

Case 2: Different Wave Packet Size (X)c= 40 m/s
LH= 30 km
τ= 12.5 min
X=25, 50, 100, 200km
Ts=120 min
Xs=400 km
dT=1 minΤs_min= 3dT
θ = 225 deg(0=eastward)• The bigger the wave packet
size (X), the sharper thespectrum• Total power increases ~4x
when X increases 2x.

Case 3: Different Wave Packet Duration (T)c= 40 m/s
LH= 30 km
τ= 12.5 min
X= 100 kmT=30, 60, 120 min
Xs=400 km
dT=1 minΤs_min= 3dT
θ = 225 deg(0=eastward)• The longer the duration ofwave packet, the spectrum
becomes sharper.
• Total power increases by ~2x
when T increases by 2x.

Aliasing Effect/Mirror artifact

c= 40 m/s
• When the target frequency is lower
LH= 30 km
than but close to the Nyquist
τ= 12.5 min
frequency, the aliasing effect still
dt=1, 3, 5, 7 min
T=120 minappeared.
Xs=400 km• It is recommended to use at leastdT=1 min
Ts_min=3dT
Τs_min= 2dT
θ = 225 deg (0=eastward)

Zero Padding Size Effecta)Zpx= 512, 1024 and 2048b)Zpt= 128, 256 and 512• MacBook Pro: dual core 2.8 GHzIntel Core i7 processor with 4 MB
cache size and 16 GB memory,while the program was singlethreaded.
• No significant differences in profile.

ü Calculation time ~ proportional with the number of zero padding
Zpx= 0.5 x default à 0.5 x decrease
Zpx= 2 x default à ~4-5 x increase
Zpt= 0.5 x default à 0.5 x decrease
Zpx= 2 x default à ~2 x increase• It should be noted that bigger zero padding give finer frequency
resolution.

Noise Sensitifity Analysis
signalSignal+white noisec= 40 m/s
LH= 30 km
τ= 12.5 min
dt=1, 3, 5, 7 min
T=120 min
Xs=400 km
dT=1 minΤs_min= 2dT
θ = 225 deg (0=eastward)• PSD enhancement caused by thewhite noise can be seen around 010 m/s in the phase spectra
density and radial profile• it is recommended to exclude thespectrum between 0-20 m/s
when the data is contaminated bynoise.

Application:
Comparison of gravity wave propagation directionsat three different latitudes

Syowa station: 9 clear-nights (April-May 2013)
Shigaraki: 10 clear-nights (April-May 2011)
Tomohon: 5 clear-nights (April-May 2016)(Perwitasari et al., 2018)14

Day-to-day variation of 2D phase velocity
spectraSyowa:
q The day-to-day variation of thephase velocity spectrum at the
Syowa Station is smaller and the
propagation direction is mainly
westward.
Shigaraki:
q The day-to-day variation of thehorizontal propagation direction islarger than that at the Syowa
Station.
Tomohon:q The day to day variation of thenightly
power
spectrummagnitude is remarkable, whichindicates the intermittency ofatmospheric gravity
waves(AGWs).
15

Average day-to-day variation of the 2D
phase velocity spectrum
Syowa:
qShows that the dominant
propagation is westward
with a phase speed <50m/s.
Shigaraki:
qShows east/southeastward
dominant propagation with aphase speed of up to ~80m/s.
Tomohon:qThe day-to-day variation in
Tomohon is too strong todiscussaveragecharacteristics; however, aphase speed of up to ~100m/s and faster is observed.
16

Summary and Future WorksØ We developed a user-friendly IDL function based on Matsuda
et al.’s (2014) 3D FFT method for airglow data analysis.
Ø This function can calculate big amount of data with reasonable
execution time.
Ø Several simulations have been done to examine theperformance analysis of this new function.
Ø This new function was applied to airglow imaging obtained
from April–May of selected years at three different latitudes:
Syowa Station (69°S, 40°E), Shigaraki (35°N, 136°E), andTomohon (1°N, 122°E), run by different groups.
Ø We demonstrated the distinct difference of the GW
propagation characteristics at three different latitudes basedon the phase velocity spectra.
Ø If you are interested in using our program, please do
not hesitate to contact us.

17