diff --git a/tools/wavetool.py b/tools/wavetool.py
index a3ddb5d4..4f4ab6de 100755
--- a/tools/wavetool.py
+++ b/tools/wavetool.py
@@ -21,6 +21,9 @@
 # is wrong with firmware scheduler.
 SMART_AMP_DELAY_THRESHOLD = 8
 
+# used to avoid zero division
+eps = 0.000000000001
+
 # Module level global variable which will store command line parameters later
 cmd = None
 
@@ -112,6 +115,8 @@ def do_wave_analysis():
         analyze_wav_smart_amp(wave, fs_wav)
     if cmd.analyze == 'wov':
         analyze_wav_wov(wave, fs_wav)
+    if cmd.analyze == 'thdn':
+        analyze_wav_thdn(wave, fs_wav)
 
 # remove digital zeros in two sides
 def trim_wave(wave):
@@ -205,6 +210,17 @@ def analyze_wav_wov(wave, fs):
         sys.exit(1002)
     print('wave analysis result: PASSED')
 
+def analyze_wav_thdn(wave, fs_wav):
+    wave = normalize(wave)
+    check_rms(wave)
+    thdn = calc_thdn(wave, fs_wav, cmd.freq[0])
+    result = thdn < cmd.threshold
+    print('THD+N of specified wave is %s dB' % thdn)
+    if not np.all(result):
+        print('THD+N too high, wave analysis result: FAILED')
+        sys.exit(1003)
+    print('wave analysis result: PASSED')
+
 def calc_thdn(wave, fs, freq):
     """
     This function calculates Total Harmonic Distortion plus Noise (THD+N) of a wave.
@@ -217,14 +233,36 @@ def calc_thdn(wave, fs, freq):
 
     Returns
     ----------
-    A list of THD+N values for each channel
+    A list of THD+N values for each channel, unit: dB
     """
     # skip 20ms to avoid the pulse due to filtering
     time_skip = 0.02
     samples_skip = int(time_skip * fs)
     filtered = stdnotch(wave, freq, fs)
-    filtered_cut = filtered[samples_skip:,:]
-    return 10 * np.log10(np.mean(np.power(filtered_cut, 2), axis=0))
+    filtered_cut = filtered[samples_skip:]
+    return 10 * np.log10(np.mean(np.power(filtered_cut, 2), axis=0) + eps)
+
+def check_rms(wave):
+    dbfs = calc_rms(wave)
+    print('RMS of specified wave is %s dB' % dbfs)
+    dbfs_result = dbfs > cmd.zero_threshold
+    if not np.all(dbfs_result):
+        print("Wave level too low, try to increase volume!!!")
+        print("Valid wave level should be bigger than %f dBFS" % cmd.zero_threshold)
+        sys.exit(1003)
+
+def calc_rms(wave):
+    """
+    This function calculates Root Mean Square (RMS) value of a wave, the value is used
+    to identify wave level.
+    Parameters
+    ----------
+    wave: input wave data
+    Returns
+    ----------
+    A list of RMS values for each channel, unit: dBFS
+    """
+    return 10 * np.log10(np.mean(np.power(wave,2), axis=0) + eps)
 
 def parse_cmdline():
     parser = argparse.ArgumentParser(add_help=True, formatter_class=argparse.RawTextHelpFormatter,
@@ -242,7 +280,7 @@ def parse_cmdline():
     help='sample bits of generated wave')
     parser.add_argument('-o', '--output', type=str, help='path to store generated files', default='.')
     # wave comparison arguments
-    parser.add_argument('-a', '--analyze', type=str, choices=['smart_amp', 'wov'], help='analyze recorded wave to give case verdict')
+    parser.add_argument('-a', '--analyze', type=str, choices=['smart_amp', 'wov', 'thdn'], help='analyze recorded wave to give case verdict')
     parser.add_argument('-R', '--recorded_wave', type=str, help='path of recorded wave')
     parser.add_argument('-Z', '--zero_threshold', type=float, default=-50.3, help='zero threshold in dBFS')
     parser.add_argument('-H', '--hb_time', type=float, default=2.1, help='history buffer size')