diff --git a/.gitignore b/.gitignore
index b0ac0d0e..d3718d6b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -28,6 +28,9 @@ coverage.xml*
 egs/*/*/data
 egs/*/*/downloads
 egs/*/*/dump
+egs/*/*/dump*
+egs/*/*/wav_dump
+egs/*/*/wav_dump*
 egs/*/*/exp
 egs/*/*/conf/tuning
 
diff --git a/egs/csd/voc1/run.sh b/egs/csd/voc1/run.sh
index fe328a0d..9be72c27 100755
--- a/egs/csd/voc1/run.sh
+++ b/egs/csd/voc1/run.sh
@@ -18,7 +18,7 @@ n_jobs=8      # number of parallel jobs in feature extraction
 conf=conf/hifigan.v1.yaml
 
 # directory path setting
-download_dir=downloads # set the directory to your database
+download_dir="/data2/tyx/dataset" # set the directory to your database
 dumpdir=dump           # directory to dump features
 
 # training related setting
diff --git a/egs/ofuton_p_utagoe_db/voc1/run.sh b/egs/ofuton_p_utagoe_db/voc1/run.sh
index ad90dcbb..f289ff82 100755
--- a/egs/ofuton_p_utagoe_db/voc1/run.sh
+++ b/egs/ofuton_p_utagoe_db/voc1/run.sh
@@ -17,7 +17,7 @@ n_jobs=8      # number of parallel jobs in feature extraction
 conf=conf/hifigan.v1.yaml
 
 # directory path setting
-download_dir=downloads # set the directory to your database
+download_dir=/data3/tyx/dataset # set the directory to your database
 dumpdir=dump           # directory to dump features
 
 # training related setting
@@ -33,7 +33,7 @@ checkpoint="" # checkpoint path to be used for decoding
 # shellcheck disable=SC1091
 . utils/parse_options.sh || exit 1;
 
-train_set="train_nodev" # name of training data directory
+train_set="train" # name of training data directory
 dev_set="dev"           # name of development data direcotry
 eval_set="eval"         # name of evaluation data direcotry
 
diff --git a/egs/opencpop/token_voc1/cmd.sh b/egs/opencpop/token_voc1/cmd.sh
new file mode 100644
index 00000000..19f34210
--- /dev/null
+++ b/egs/opencpop/token_voc1/cmd.sh
@@ -0,0 +1,91 @@
+# ====== About run.pl, queue.pl, slurm.pl, and ssh.pl ======
+# Usage: <cmd>.pl [options] JOB=1:<nj> <log> <command...>
+# e.g.
+#   run.pl --mem 4G JOB=1:10 echo.JOB.log echo JOB
+#
+# Options:
+#   --time <time>: Limit the maximum time to execute.
+#   --mem <mem>: Limit the maximum memory usage.
+#   -–max-jobs-run <njob>: Limit the number parallel jobs. This is ignored for non-array jobs.
+#   --num-threads <ngpu>: Specify the number of CPU core.
+#   --gpu <ngpu>: Specify the number of GPU devices.
+#   --config: Change the configuration file from default.
+#
+# "JOB=1:10" is used for "array jobs" and it can control the number of parallel jobs.
+# The left string of "=", i.e. "JOB", is replaced by <N>(Nth job) in the command and the log file name,
+# e.g. "echo JOB" is changed to "echo 3" for the 3rd job and "echo 8" for 8th job respectively.
+# Note that the number must start with a positive number, so you can't use "JOB=0:10" for example.
+#
+# run.pl, queue.pl, slurm.pl, and ssh.pl have unified interface, not depending on its backend.
+# These options are mapping to specific options for each backend and
+# it is configured by "conf/queue.conf" and "conf/slurm.conf" by default.
+# If jobs failed, your configuration might be wrong for your environment.
+#
+#
+# The official documentaion for run.pl, queue.pl, slurm.pl, and ssh.pl:
+#   "Parallelization in Kaldi": http://kaldi-asr.org/doc/queue.html
+# =========================================================~
+
+
+# Select the backend used by run.sh from "local", "stdout", "sge", "slurm", or "ssh"
+cmd_backend="local"
+
+# Local machine, without any Job scheduling system
+if [ "${cmd_backend}" = local ]; then
+
+    # The other usage
+    export train_cmd="utils/run.pl"
+    # Used for "*_train.py": "--gpu" is appended optionally by run.sh
+    export cuda_cmd="utils/run.pl"
+    # Used for "*_recog.py"
+    export decode_cmd="utils/run.pl"
+
+# Local machine, without any Job scheduling system
+elif [ "${cmd_backend}" = stdout ]; then
+
+    # The other usage
+    export train_cmd="utils/stdout.pl"
+    # Used for "*_train.py": "--gpu" is appended optionally by run.sh
+    export cuda_cmd="utils/stdout.pl"
+    # Used for "*_recog.py"
+    export decode_cmd="utils/stdout.pl"
+
+# "qsub" (SGE, Torque, PBS, etc.)
+elif [ "${cmd_backend}" = sge ]; then
+    # The default setting is written in conf/queue.conf.
+    # You must change "-q g.q" for the "queue" for your environment.
+    # To know the "queue" names, type "qhost -q"
+    # Note that to use "--gpu *", you have to setup "complex_value" for the system scheduler.
+
+    export train_cmd="utils/queue.pl"
+    export cuda_cmd="utils/queue.pl"
+    export decode_cmd="utils/queue.pl"
+
+# "sbatch" (Slurm)
+elif [ "${cmd_backend}" = slurm ]; then
+    # The default setting is written in conf/slurm.conf.
+    # You must change "-p cpu" and "-p gpu" for the "partion" for your environment.
+    # To know the "partion" names, type "sinfo".
+    # You can use "--gpu * " by defualt for slurm and it is interpreted as "--gres gpu:*"
+    # The devices are allocated exclusively using "${CUDA_VISIBLE_DEVICES}".
+
+    export train_cmd="utils/slurm.pl"
+    export cuda_cmd="utils/slurm.pl"
+    export decode_cmd="utils/slurm.pl"
+
+elif [ "${cmd_backend}" = ssh ]; then
+    # You have to create ".queue/machines" to specify the host to execute jobs.
+    # e.g. .queue/machines
+    #   host1
+    #   host2
+    #   host3
+    # Assuming you can login them without any password, i.e. You have to set ssh keys.
+
+    export train_cmd="utils/ssh.pl"
+    export cuda_cmd="utils/ssh.pl"
+    export decode_cmd="utils/ssh.pl"
+
+else
+    echo "$0: Error: Unknown cmd_backend=${cmd_backend}" 1>&2
+    return 1
+fi
diff --git a/egs/opencpop/token_voc1/conf/hifigan_token_16k_duration.v1.yaml b/egs/opencpop/token_voc1/conf/hifigan_token_16k_duration.v1.yaml
new file mode 100644
index 00000000..68e9e4a9
--- /dev/null
+++ b/egs/opencpop/token_voc1/conf/hifigan_token_16k_duration.v1.yaml
@@ -0,0 +1,185 @@
+# This configuration is based on HiFiGAN V1, derived
+# from official repository (https://github.com/jik876/hifi-gan).
+
+# for 16k pre-trained model token, like HuBERT, WavLM, XSL-R...
+
+###########################################################
+#                FEATURE EXTRACTION SETTING               #
+###########################################################
+sampling_rate: 16000     # Sampling rate.
+fft_size: null           # FFT size.
+hop_size: 320            # Hop size.
+win_length: null         # Window length.
+                         # If set to null, it will be the same as fft_size.
+window: null             # Window function.
+num_mels: 1              # Number of mel basis.
+fmin: null               # Minimum freq in mel basis calculation.
+fmax: null               # Maximum frequency in mel basis calculation.
+global_gain_scale: 1.0   # Will be multiplied to all of waveform.
+trim_silence: false      # Whether to trim the start and end of silence.
+trim_threshold_in_db: 20 # Need to tune carefully if the recording is not good.
+trim_frame_size: 1024    # Frame size in trimming.
+trim_hop_size: 256       # Hop size in trimming.
+format: "hdf5"           # Feature file format. "npy" or "hdf5" is supported.
+
+###########################################################
+#         GENERATOR NETWORK ARCHITECTURE SETTING          #
+###########################################################
+generator_type: DiscreteSymbolDurationGenerator
+generator_params:
+    in_channels: 1024                      # Number of input channels.
+    out_channels: 1                       # Number of output channels.
+    channels: 512                         # Number of initial channels.
+    num_embs: 1024
+    kernel_size: 7                        # Kernel size of initial and final conv layers.
+    upsample_scales: [10, 8, 2, 2]        # Upsampling scales.
+    upsample_kernel_sizes: [20, 16, 4, 4] # Kernel size for upsampling layers.
+    resblock_kernel_sizes: [3, 7, 11]     # Kernel size for residual blocks.
+    resblock_dilations:                   # Dilations for residual blocks.
+        - [1, 3, 5]
+        - [1, 3, 5]
+        - [1, 3, 5]
+    use_additional_convs: true            # Whether to use additional conv layer in residual blocks.
+    bias: true                            # Whether to use bias parameter in conv.
+    nonlinear_activation: "LeakyReLU"     # Nonlinear activation type.
+    nonlinear_activation_params:          # Nonlinear activation paramters.
+        negative_slope: 0.1
+    use_weight_norm: true                 # Whether to apply weight normalization.
+    duration_layers: 2                    # Duration predictor layers
+    duration_chans: 384                   # Duration predictor conv channels
+    duration_kernel_size: 3               # Duration predictor kernel size
+    duration_offset: 1.0                  # Duration predictor offset 
+    duration_dropout_rate: 0.5            # Duration predictor dropout
+    num_spk_embs: 0                       # Do not consider speaker embedding for single spk
+
+###########################################################
+#       DISCRIMINATOR NETWORK ARCHITECTURE SETTING        #
+###########################################################
+discriminator_type: HiFiGANMultiScaleMultiPeriodDiscriminator
+discriminator_params:
+    scales: 3                              # Number of multi-scale discriminator.
+    scale_downsample_pooling: "AvgPool1d"  # Pooling operation for scale discriminator.
+    scale_downsample_pooling_params:
+        kernel_size: 4                     # Pooling kernel size.
+        stride: 2                          # Pooling stride.
+        padding: 2                         # Padding size.
+    scale_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [15, 41, 5, 3]       # List of kernel sizes.
+        channels: 128                      # Initial number of channels.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        max_groups: 16                     # Maximum number of groups in downsampling conv layers.
+        bias: true
+        downsample_scales: [4, 4, 4, 4, 1] # Downsampling scales.
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:
+            negative_slope: 0.1
+    follow_official_norm: true             # Whether to follow the official norm setting.
+    periods: [2, 3, 5, 7, 11]              # List of period for multi-period discriminator.
+    period_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [5, 3]               # List of kernel sizes.
+        channels: 32                       # Initial number of channels.
+        downsample_scales: [3, 3, 3, 3, 1] # Downsampling scales.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        bias: true                         # Whether to use bias parameter in conv layer."
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:       # Nonlinear activation paramters.
+            negative_slope: 0.1
+        use_weight_norm: true              # Whether to apply weight normalization.
+        use_spectral_norm: false           # Whether to apply spectral normalization.
+
+###########################################################
+#                   STFT LOSS SETTING                     #
+###########################################################
+use_stft_loss: false # Whether to use multi-resolution STFT loss.
+use_duration_loss: true # Whether to use duration prediction loss (need for duration predictor)
+duration_loss_params:
+    offset: 1.0
+    reduction: mean
+use_mel_loss: true   # Whether to use Mel-spectrogram loss.
+mel_loss_params:     # Mel-spectrogram loss parameters.
+    fs: 16000
+    fft_size: 1024
+    hop_size: 256
+    win_length: null
+    window: "hann"
+    num_mels: 80
+    fmin: 0
+    fmax: 8000
+    log_base: null   # Log base. If set to null, use natural logarithm.
+generator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+discriminator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+use_feat_match_loss: true
+feat_match_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+    average_by_layers: false         # Whether to average loss by #layers in each discriminator.
+    include_final_outputs: true      # Whether to include final outputs in feat match loss calculation.
+
+###########################################################
+#               ADVERSARIAL LOSS SETTING                  #
+###########################################################
+lambda_aux: 45.0       # Loss balancing coefficient for STFT loss.
+lambda_adv: 1.0        # Loss balancing coefficient for adversarial loss.
+lambda_feat_match: 2.0 # Loss balancing coefficient for feat match loss..
+
+###########################################################
+#                  DATA LOADER SETTING                    #
+###########################################################
+batch_size: 16              # Batch size.
+batch_max_steps: 10240      # Length of each audio in batch. Make sure dividable by hop_size.
+pin_memory: true            # Whether to pin memory in Pytorch DataLoader.
+num_workers: 0              # Number of workers in Pytorch DataLoader.
+remove_short_samples: false # Whether to remove samples the length of which are less than batch_max_steps.
+allow_cache: true           # Whether to allow cache in dataset. If true, it requires cpu memory.
+
+###########################################################
+#             OPTIMIZER & SCHEDULER SETTING               #
+###########################################################
+generator_optimizer_type: Adam
+generator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+generator_scheduler_type: MultiStepLR
+generator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+generator_grad_norm: -1
+discriminator_optimizer_type: Adam
+discriminator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+discriminator_scheduler_type: MultiStepLR
+discriminator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+discriminator_grad_norm: -1
+
+###########################################################
+#                    INTERVAL SETTING                     #
+###########################################################
+generator_train_start_steps: 1     # Number of steps to start to train discriminator.
+discriminator_train_start_steps: 0 # Number of steps to start to train discriminator.
+train_max_steps: 2500000           # Number of training steps.
+save_interval_steps: 50000         # Interval steps to save checkpoint.
+eval_interval_steps: 1000          # Interval steps to evaluate the network.
+log_interval_steps: 100            # Interval steps to record the training log.
+
+###########################################################
+#                     OTHER SETTING                       #
+###########################################################
+num_save_intermediate_results: 4  # Number of results to be saved as intermediate results.
diff --git a/egs/opencpop/token_voc1/conf/hifigan_token_16k_nodp.v1.yaml b/egs/opencpop/token_voc1/conf/hifigan_token_16k_nodp.v1.yaml
new file mode 100644
index 00000000..59d510b1
--- /dev/null
+++ b/egs/opencpop/token_voc1/conf/hifigan_token_16k_nodp.v1.yaml
@@ -0,0 +1,180 @@
+# This configuration is based on HiFiGAN V1, derived
+# from official repository (https://github.com/jik876/hifi-gan).
+
+# for 16k pre-trained model token, like HuBERT, WavLM, XSL-R...
+
+###########################################################
+#                FEATURE EXTRACTION SETTING               #
+###########################################################
+sampling_rate: 16000     # Sampling rate.
+fft_size: null           # FFT size.
+hop_size: 320            # Hop size.
+win_length: null         # Window length.
+                         # If set to null, it will be the same as fft_size.
+window: null             # Window function.
+num_mels: 1              # Number of mel basis.
+fmin: null               # Minimum freq in mel basis calculation.
+fmax: null               # Maximum frequency in mel basis calculation.
+global_gain_scale: 1.0   # Will be multiplied to all of waveform.
+trim_silence: false      # Whether to trim the start and end of silence.
+trim_threshold_in_db: 20 # Need to tune carefully if the recording is not good.
+trim_frame_size: 1024    # Frame size in trimming.
+trim_hop_size: 256       # Hop size in trimming.
+format: "hdf5"           # Feature file format. "npy" or "hdf5" is supported.
+
+###########################################################
+#         GENERATOR NETWORK ARCHITECTURE SETTING          #
+###########################################################
+generator_type: DiscreteSymbolHiFiGANGenerator
+generator_params:
+    in_channels: 768                      # Number of input channels.
+    out_channels: 1                       # Number of output channels.
+    channels: 512                         # Number of initial channels.
+    num_embs: 1025
+    kernel_size: 7                        # Kernel size of initial and final conv layers.
+    upsample_scales: [10, 8, 2, 2]        # Upsampling scales.
+    upsample_kernel_sizes: [20, 16, 4, 4] # Kernel size for upsampling layers.
+    resblock_kernel_sizes: [3, 7, 11]     # Kernel size for residual blocks.
+    resblock_dilations:                   # Dilations for residual blocks.
+        - [1, 3, 5]
+        - [1, 3, 5]
+        - [1, 3, 5]
+    use_additional_convs: true            # Whether to use additional conv layer in residual blocks.
+    bias: true                            # Whether to use bias parameter in conv.
+    nonlinear_activation: "LeakyReLU"     # Nonlinear activation type.
+    nonlinear_activation_params:          # Nonlinear activation paramters.
+        negative_slope: 0.1
+    use_weight_norm: true                 # Whether to apply weight normalization.
+    num_spk_embs: 0                       # Do not consider speaker embedding for single spk
+
+###########################################################
+#       DISCRIMINATOR NETWORK ARCHITECTURE SETTING        #
+###########################################################
+discriminator_type: HiFiGANMultiScaleMultiPeriodDiscriminator
+discriminator_params:
+    scales: 3                              # Number of multi-scale discriminator.
+    scale_downsample_pooling: "AvgPool1d"  # Pooling operation for scale discriminator.
+    scale_downsample_pooling_params:
+        kernel_size: 4                     # Pooling kernel size.
+        stride: 2                          # Pooling stride.
+        padding: 2                         # Padding size.
+    scale_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [15, 41, 5, 3]       # List of kernel sizes.
+        channels: 128                      # Initial number of channels.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        max_groups: 16                     # Maximum number of groups in downsampling conv layers.
+        bias: true
+        downsample_scales: [4, 4, 4, 4, 1] # Downsampling scales.
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:
+            negative_slope: 0.1
+    follow_official_norm: true             # Whether to follow the official norm setting.
+    periods: [2, 3, 5, 7, 11]              # List of period for multi-period discriminator.
+    period_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [5, 3]               # List of kernel sizes.
+        channels: 32                       # Initial number of channels.
+        downsample_scales: [3, 3, 3, 3, 1] # Downsampling scales.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        bias: true                         # Whether to use bias parameter in conv layer."
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:       # Nonlinear activation paramters.
+            negative_slope: 0.1
+        use_weight_norm: true              # Whether to apply weight normalization.
+        use_spectral_norm: false           # Whether to apply spectral normalization.
+
+###########################################################
+#                   STFT LOSS SETTING                     #
+###########################################################
+use_stft_loss: false # Whether to use multi-resolution STFT loss.
+use_duration_loss: true # Whether to use duration prediction loss (need for duration predictor)
+duration_loss_params:
+    offset: 1.0
+    reduction: mean
+use_mel_loss: true   # Whether to use Mel-spectrogram loss.
+mel_loss_params:     # Mel-spectrogram loss parameters.
+    fs: 16000
+    fft_size: 1024
+    hop_size: 256
+    win_length: null
+    window: "hann"
+    num_mels: 80
+    fmin: 0
+    fmax: 8000
+    log_base: null   # Log base. If set to null, use natural logarithm.
+generator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+discriminator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+use_feat_match_loss: true
+feat_match_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+    average_by_layers: false         # Whether to average loss by #layers in each discriminator.
+    include_final_outputs: true      # Whether to include final outputs in feat match loss calculation.
+
+###########################################################
+#               ADVERSARIAL LOSS SETTING                  #
+###########################################################
+lambda_aux: 45.0       # Loss balancing coefficient for STFT loss.
+lambda_adv: 1.0        # Loss balancing coefficient for adversarial loss.
+lambda_feat_match: 2.0 # Loss balancing coefficient for feat match loss..
+
+###########################################################
+#                  DATA LOADER SETTING                    #
+###########################################################
+batch_size: 16              # Batch size.
+batch_max_steps: 10240      # Length of each audio in batch. Make sure dividable by hop_size.
+pin_memory: true            # Whether to pin memory in Pytorch DataLoader.
+num_workers: 0              # Number of workers in Pytorch DataLoader.
+remove_short_samples: false # Whether to remove samples the length of which are less than batch_max_steps.
+allow_cache: true           # Whether to allow cache in dataset. If true, it requires cpu memory.
+
+###########################################################
+#             OPTIMIZER & SCHEDULER SETTING               #
+###########################################################
+generator_optimizer_type: Adam
+generator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+generator_scheduler_type: MultiStepLR
+generator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+generator_grad_norm: -1
+discriminator_optimizer_type: Adam
+discriminator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+discriminator_scheduler_type: MultiStepLR
+discriminator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+discriminator_grad_norm: -1
+
+###########################################################
+#                    INTERVAL SETTING                     #
+###########################################################
+generator_train_start_steps: 1     # Number of steps to start to train discriminator.
+discriminator_train_start_steps: 0 # Number of steps to start to train discriminator.
+train_max_steps: 400000           # Number of training steps.
+save_interval_steps: 50000         # Interval steps to save checkpoint.
+eval_interval_steps: 10000          # Interval steps to evaluate the network.
+log_interval_steps: 100            # Interval steps to record the training log.
+
+###########################################################
+#                     OTHER SETTING                       #
+###########################################################
+num_save_intermediate_results: 4  # Number of results to be saved as intermediate results.
diff --git a/egs/opencpop/token_voc1/conf/hifigan_token_16k_nodp_f0.v1.yaml b/egs/opencpop/token_voc1/conf/hifigan_token_16k_nodp_f0.v1.yaml
new file mode 100644
index 00000000..d5d2c5e3
--- /dev/null
+++ b/egs/opencpop/token_voc1/conf/hifigan_token_16k_nodp_f0.v1.yaml
@@ -0,0 +1,187 @@
+# This configuration is based on HiFiGAN V1, derived
+# from official repository (https://github.com/jik876/hifi-gan).
+
+# for 16k pre-trained model token, like HuBERT, WavLM, XSL-R...
+
+###########################################################
+#                FEATURE EXTRACTION SETTING               #
+###########################################################
+sampling_rate: 16000     # Sampling rate.
+fft_size: null           # FFT size.
+hop_size: 320            # Hop size.
+win_length: null         # Window length.
+                         # If set to null, it will be the same as fft_size.
+window: null             # Window function.
+num_mels: 1              # Number of mel basis.
+fmin: null               # Minimum freq in mel basis calculation.
+fmax: null               # Maximum frequency in mel basis calculation.
+global_gain_scale: 1.0   # Will be multiplied to all of waveform.
+trim_silence: false      # Whether to trim the start and end of silence.
+trim_threshold_in_db: 20 # Need to tune carefully if the recording is not good.
+trim_frame_size: 1024    # Frame size in trimming.
+trim_hop_size: 256       # Hop size in trimming.
+format: "hdf5"           # Feature file format. "npy" or "hdf5" is supported.
+
+###########################################################
+#         GENERATOR NETWORK ARCHITECTURE SETTING          #
+###########################################################
+generator_type: DiscreteSymbolF0Generator
+generator_params:
+    in_channels: 768                      # Number of input channels.
+    out_channels: 1                       # Number of output channels.
+    channels: 512                         # Number of initial channels.
+    linear_channel: 256                   # Number of Linear channels(f0).
+    num_embs: 1025
+    kernel_size: 7                        # Kernel size of initial and final conv layers.
+    upsample_scales: [10, 8, 2, 2]        # Upsampling scales.
+    upsample_kernel_sizes: [20, 16, 4, 4] # Kernel size for upsampling layers.
+    resblock_kernel_sizes: [3, 7, 11]     # Kernel size for residual blocks.
+    resblock_dilations:                   # Dilations for residual blocks.
+        - [1, 3, 5]
+        - [1, 3, 5]
+        - [1, 3, 5]
+    use_additional_convs: true            # Whether to use additional conv layer in residual blocks.
+    bias: true                            # Whether to use bias parameter in conv.
+    nonlinear_activation: "LeakyReLU"     # Nonlinear activation type.
+    nonlinear_activation_params:          # Nonlinear activation paramters.
+        negative_slope: 0.1
+    use_weight_norm: true                 # Whether to apply weight normalization.
+    num_spk_embs: 0                       # Do not consider speaker embedding for single spk
+    
+    use_embedding_feats: false            # Whether to use embedding feature of models
+    use_weight_sum: false                  # Whether to apply weight sum
+    layer_num: 12                         # Number of layers used
+    use_fix_weight: false                 # Whether to fix parameter of weight sum (for Residual cluster layer)
+    use_f0: true                          # Whether to use f0
+
+###########################################################
+#       DISCRIMINATOR NETWORK ARCHITECTURE SETTING        #
+###########################################################
+discriminator_type: HiFiGANMultiScaleMultiPeriodDiscriminator
+discriminator_params:
+    scales: 3                              # Number of multi-scale discriminator.
+    scale_downsample_pooling: "AvgPool1d"  # Pooling operation for scale discriminator.
+    scale_downsample_pooling_params:
+        kernel_size: 4                     # Pooling kernel size.
+        stride: 2                          # Pooling stride.
+        padding: 2                         # Padding size.
+    scale_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [15, 41, 5, 3]       # List of kernel sizes.
+        channels: 128                      # Initial number of channels.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        max_groups: 16                     # Maximum number of groups in downsampling conv layers.
+        bias: true
+        downsample_scales: [4, 4, 4, 4, 1] # Downsampling scales.
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:
+            negative_slope: 0.1
+    follow_official_norm: true             # Whether to follow the official norm setting.
+    periods: [2, 3, 5, 7, 11]              # List of period for multi-period discriminator.
+    period_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [5, 3]               # List of kernel sizes.
+        channels: 32                       # Initial number of channels.
+        downsample_scales: [3, 3, 3, 3, 1] # Downsampling scales.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        bias: true                         # Whether to use bias parameter in conv layer."
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:       # Nonlinear activation paramters.
+            negative_slope: 0.1
+        use_weight_norm: true              # Whether to apply weight normalization.
+        use_spectral_norm: false           # Whether to apply spectral normalization.
+
+###########################################################
+#                   STFT LOSS SETTING                     #
+###########################################################
+use_stft_loss: false # Whether to use multi-resolution STFT loss.
+use_duration_loss: true # Whether to use duration prediction loss (need for duration predictor)
+duration_loss_params:
+    offset: 1.0
+    reduction: mean
+use_mel_loss: true   # Whether to use Mel-spectrogram loss.
+mel_loss_params:     # Mel-spectrogram loss parameters.
+    fs: 16000
+    fft_size: 1024
+    hop_size: 256
+    win_length: null
+    window: "hann"
+    num_mels: 80
+    fmin: 0
+    fmax: 8000
+    log_base: null   # Log base. If set to null, use natural logarithm.
+generator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+discriminator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+use_feat_match_loss: true
+feat_match_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+    average_by_layers: false         # Whether to average loss by #layers in each discriminator.
+    include_final_outputs: true      # Whether to include final outputs in feat match loss calculation.
+
+###########################################################
+#               ADVERSARIAL LOSS SETTING                  #
+###########################################################
+lambda_aux: 45.0       # Loss balancing coefficient for STFT loss.
+lambda_adv: 1.0        # Loss balancing coefficient for adversarial loss.
+lambda_feat_match: 2.0 # Loss balancing coefficient for feat match loss..
+
+###########################################################
+#                  DATA LOADER SETTING                    #
+###########################################################
+batch_size: 16              # Batch size.
+batch_max_steps: 10240      # Length of each audio in batch. Make sure dividable by hop_size.
+pin_memory: true            # Whether to pin memory in Pytorch DataLoader.
+num_workers: 0              # Number of workers in Pytorch DataLoader.
+remove_short_samples: false # Whether to remove samples the length of which are less than batch_max_steps.
+allow_cache: true           # Whether to allow cache in dataset. If true, it requires cpu memory.
+
+###########################################################
+#             OPTIMIZER & SCHEDULER SETTING               #
+###########################################################
+generator_optimizer_type: Adam
+generator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+generator_scheduler_type: MultiStepLR
+generator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+generator_grad_norm: -1
+discriminator_optimizer_type: Adam
+discriminator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+discriminator_scheduler_type: MultiStepLR
+discriminator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+discriminator_grad_norm: -1
+
+###########################################################
+#                    INTERVAL SETTING                     #
+###########################################################
+generator_train_start_steps: 1     # Number of steps to start to train discriminator.
+discriminator_train_start_steps: 0 # Number of steps to start to train discriminator.
+train_max_steps: 250000           # Number of training steps.
+save_interval_steps: 50000         # Interval steps to save checkpoint.
+eval_interval_steps: 10000          # Interval steps to evaluate the network.
+log_interval_steps: 100            # Interval steps to record the training log.
+
+###########################################################
+#                     OTHER SETTING                       #
+###########################################################
+num_save_intermediate_results: 4  # Number of results to be saved as intermediate results.
diff --git a/egs/opencpop/token_voc1/conf/hifigan_token_24k_nodp_f0.v1.yaml b/egs/opencpop/token_voc1/conf/hifigan_token_24k_nodp_f0.v1.yaml
new file mode 100644
index 00000000..99e5da55
--- /dev/null
+++ b/egs/opencpop/token_voc1/conf/hifigan_token_24k_nodp_f0.v1.yaml
@@ -0,0 +1,187 @@
+# This configuration is based on HiFiGAN V1, derived
+# from official repository (https://github.com/jik876/hifi-gan).
+
+# for 24k pre-trained model token, like MERT...
+
+###########################################################
+#                FEATURE EXTRACTION SETTING               #
+###########################################################
+sampling_rate: 24000     # Sampling rate.
+fft_size: null           # FFT size.
+hop_size: 320            # Hop size.
+win_length: null         # Window length.
+                         # If set to null, it will be the same as fft_size.
+window: null             # Window function.
+num_mels: 1              # Number of mel basis.
+fmin: null               # Minimum freq in mel basis calculation.
+fmax: null               # Maximum frequency in mel basis calculation.
+global_gain_scale: 1.0   # Will be multiplied to all of waveform.
+trim_silence: false      # Whether to trim the start and end of silence.
+trim_threshold_in_db: 20 # Need to tune carefully if the recording is not good.
+trim_frame_size: 1024    # Frame size in trimming.
+trim_hop_size: 256       # Hop size in trimming.
+format: "hdf5"           # Feature file format. "npy" or "hdf5" is supported.
+
+###########################################################
+#         GENERATOR NETWORK ARCHITECTURE SETTING          #
+###########################################################
+generator_type: DiscreteSymbolF0Generator
+generator_params:
+    in_channels: 768                      # Number of input channels.
+    out_channels: 1                       # Number of output channels.
+    channels: 512                         # Number of initial channels.
+    linear_channel: 256                   # Number of Linear channels(f0).
+    num_embs: 1025
+    kernel_size: 7                        # Kernel size of initial and final conv layers.
+    upsample_scales: [10, 8, 2, 2]        # Upsampling scales.
+    upsample_kernel_sizes: [20, 16, 4, 4] # Kernel size for upsampling layers.
+    resblock_kernel_sizes: [3, 7, 11]     # Kernel size for residual blocks.
+    resblock_dilations:                   # Dilations for residual blocks.
+        - [1, 3, 5]
+        - [1, 3, 5]
+        - [1, 3, 5]
+    use_additional_convs: true            # Whether to use additional conv layer in residual blocks.
+    bias: true                            # Whether to use bias parameter in conv.
+    nonlinear_activation: "LeakyReLU"     # Nonlinear activation type.
+    nonlinear_activation_params:          # Nonlinear activation paramters.
+        negative_slope: 0.1
+    use_weight_norm: true                 # Whether to apply weight normalization.
+    num_spk_embs: 0                       # Do not consider speaker embedding for single spk
+    
+    use_embedding_feats: false            # Whether to use embedding feature of models
+    use_weight_sum: false                  # Whether to apply weight sum
+    layer_num: 12                         # Number of layers used
+    use_fix_weight: false                 # Whether to fix parameter of weight sum (for Residual cluster layer)
+    use_f0: true                          # Whether to use f0
+
+###########################################################
+#       DISCRIMINATOR NETWORK ARCHITECTURE SETTING        #
+###########################################################
+discriminator_type: HiFiGANMultiScaleMultiPeriodDiscriminator
+discriminator_params:
+    scales: 3                              # Number of multi-scale discriminator.
+    scale_downsample_pooling: "AvgPool1d"  # Pooling operation for scale discriminator.
+    scale_downsample_pooling_params:
+        kernel_size: 4                     # Pooling kernel size.
+        stride: 2                          # Pooling stride.
+        padding: 2                         # Padding size.
+    scale_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [15, 41, 5, 3]       # List of kernel sizes.
+        channels: 128                      # Initial number of channels.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        max_groups: 16                     # Maximum number of groups in downsampling conv layers.
+        bias: true
+        downsample_scales: [4, 4, 4, 4, 1] # Downsampling scales.
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:
+            negative_slope: 0.1
+    follow_official_norm: true             # Whether to follow the official norm setting.
+    periods: [2, 3, 5, 7, 11]              # List of period for multi-period discriminator.
+    period_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [5, 3]               # List of kernel sizes.
+        channels: 32                       # Initial number of channels.
+        downsample_scales: [3, 3, 3, 3, 1] # Downsampling scales.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        bias: true                         # Whether to use bias parameter in conv layer."
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:       # Nonlinear activation paramters.
+            negative_slope: 0.1
+        use_weight_norm: true              # Whether to apply weight normalization.
+        use_spectral_norm: false           # Whether to apply spectral normalization.
+
+###########################################################
+#                   STFT LOSS SETTING                     #
+###########################################################
+use_stft_loss: false # Whether to use multi-resolution STFT loss.
+use_duration_loss: true # Whether to use duration prediction loss (need for duration predictor)
+duration_loss_params:
+    offset: 1.0
+    reduction: mean
+use_mel_loss: true   # Whether to use Mel-spectrogram loss.
+mel_loss_params:     # Mel-spectrogram loss parameters.
+    fs: 24000
+    fft_size: 2048
+    hop_size: 300
+    win_length: 1200
+    window: "hann"
+    num_mels: 80
+    fmin: 0
+    fmax: 8000
+    log_base: null   # Log base. If set to null, use natural logarithm.
+generator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+discriminator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+use_feat_match_loss: true
+feat_match_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+    average_by_layers: false         # Whether to average loss by #layers in each discriminator.
+    include_final_outputs: true      # Whether to include final outputs in feat match loss calculation.
+
+###########################################################
+#               ADVERSARIAL LOSS SETTING                  #
+###########################################################
+lambda_aux: 45.0       # Loss balancing coefficient for STFT loss.
+lambda_adv: 1.0        # Loss balancing coefficient for adversarial loss.
+lambda_feat_match: 2.0 # Loss balancing coefficient for feat match loss..
+
+###########################################################
+#                  DATA LOADER SETTING                    #
+###########################################################
+batch_size: 16              # Batch size.
+batch_max_steps: 10240      # Length of each audio in batch. Make sure dividable by hop_size.
+pin_memory: true            # Whether to pin memory in Pytorch DataLoader.
+num_workers: 0              # Number of workers in Pytorch DataLoader.
+remove_short_samples: false # Whether to remove samples the length of which are less than batch_max_steps.
+allow_cache: true           # Whether to allow cache in dataset. If true, it requires cpu memory.
+
+###########################################################
+#             OPTIMIZER & SCHEDULER SETTING               #
+###########################################################
+generator_optimizer_type: Adam
+generator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+generator_scheduler_type: MultiStepLR
+generator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+generator_grad_norm: -1
+discriminator_optimizer_type: Adam
+discriminator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+discriminator_scheduler_type: MultiStepLR
+discriminator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+discriminator_grad_norm: -1
+
+###########################################################
+#                    INTERVAL SETTING                     #
+###########################################################
+generator_train_start_steps: 1     # Number of steps to start to train discriminator.
+discriminator_train_start_steps: 0 # Number of steps to start to train discriminator.
+train_max_steps: 250000           # Number of training steps.
+save_interval_steps: 50000         # Interval steps to save checkpoint.
+eval_interval_steps: 10000          # Interval steps to evaluate the network.
+log_interval_steps: 100            # Interval steps to record the training log.
+
+###########################################################
+#                     OTHER SETTING                       #
+###########################################################
+num_save_intermediate_results: 4  # Number of results to be saved as intermediate results.
diff --git a/egs/opencpop/token_voc1/conf/hifigan_token_48k_nodp.v1.yaml b/egs/opencpop/token_voc1/conf/hifigan_token_48k_nodp.v1.yaml
new file mode 100644
index 00000000..ae4ea68d
--- /dev/null
+++ b/egs/opencpop/token_voc1/conf/hifigan_token_48k_nodp.v1.yaml
@@ -0,0 +1,180 @@
+# This configuration is based on HiFiGAN V1, derived
+# from official repository (https://github.com/jik876/hifi-gan).
+
+# for 16k pre-trained model token, like HuBERT, WavLM, XSL-R...
+
+###########################################################
+#                FEATURE EXTRACTION SETTING               #
+###########################################################
+sampling_rate: 48000     # Sampling rate.
+fft_size: null           # FFT size.
+hop_size: 320            # Hop size.
+win_length: null         # Window length.
+                         # If set to null, it will be the same as fft_size.
+window: null             # Window function.
+num_mels: 1              # Number of mel basis.
+fmin: null               # Minimum freq in mel basis calculation.
+fmax: null               # Maximum frequency in mel basis calculation.
+global_gain_scale: 1.0   # Will be multiplied to all of waveform.
+trim_silence: false      # Whether to trim the start and end of silence.
+trim_threshold_in_db: 20 # Need to tune carefully if the recording is not good.
+trim_frame_size: 1024    # Frame size in trimming.
+trim_hop_size: 256       # Hop size in trimming.
+format: "hdf5"           # Feature file format. "npy" or "hdf5" is supported.
+
+###########################################################
+#         GENERATOR NETWORK ARCHITECTURE SETTING          #
+###########################################################
+generator_type: DiscreteSymbolHiFiGANGenerator
+generator_params:
+    in_channels: 768                      # Number of input channels.
+    out_channels: 1                       # Number of output channels.
+    channels: 512                         # Number of initial channels.
+    num_embs: 1025
+    kernel_size: 7                        # Kernel size of initial and final conv layers.
+    upsample_scales: [10, 8, 2, 2]        # Upsampling scales.
+    upsample_kernel_sizes: [20, 16, 4, 4] # Kernel size for upsampling layers.
+    resblock_kernel_sizes: [3, 7, 11]     # Kernel size for residual blocks.
+    resblock_dilations:                   # Dilations for residual blocks.
+        - [1, 3, 5]
+        - [1, 3, 5]
+        - [1, 3, 5]
+    use_additional_convs: true            # Whether to use additional conv layer in residual blocks.
+    bias: true                            # Whether to use bias parameter in conv.
+    nonlinear_activation: "LeakyReLU"     # Nonlinear activation type.
+    nonlinear_activation_params:          # Nonlinear activation paramters.
+        negative_slope: 0.1
+    use_weight_norm: true                 # Whether to apply weight normalization.
+    num_spk_embs: 0                       # Do not consider speaker embedding for single spk
+
+###########################################################
+#       DISCRIMINATOR NETWORK ARCHITECTURE SETTING        #
+###########################################################
+discriminator_type: HiFiGANMultiScaleMultiPeriodDiscriminator
+discriminator_params:
+    scales: 3                              # Number of multi-scale discriminator.
+    scale_downsample_pooling: "AvgPool1d"  # Pooling operation for scale discriminator.
+    scale_downsample_pooling_params:
+        kernel_size: 4                     # Pooling kernel size.
+        stride: 2                          # Pooling stride.
+        padding: 2                         # Padding size.
+    scale_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [15, 41, 5, 3]       # List of kernel sizes.
+        channels: 128                      # Initial number of channels.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        max_groups: 16                     # Maximum number of groups in downsampling conv layers.
+        bias: true
+        downsample_scales: [4, 4, 4, 4, 1] # Downsampling scales.
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:
+            negative_slope: 0.1
+    follow_official_norm: true             # Whether to follow the official norm setting.
+    periods: [2, 3, 5, 7, 11]              # List of period for multi-period discriminator.
+    period_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [5, 3]               # List of kernel sizes.
+        channels: 32                       # Initial number of channels.
+        downsample_scales: [3, 3, 3, 3, 1] # Downsampling scales.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        bias: true                         # Whether to use bias parameter in conv layer."
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:       # Nonlinear activation paramters.
+            negative_slope: 0.1
+        use_weight_norm: true              # Whether to apply weight normalization.
+        use_spectral_norm: false           # Whether to apply spectral normalization.
+
+###########################################################
+#                   STFT LOSS SETTING                     #
+###########################################################
+use_stft_loss: false # Whether to use multi-resolution STFT loss.
+use_duration_loss: true # Whether to use duration prediction loss (need for duration predictor)
+duration_loss_params:
+    offset: 1.0
+    reduction: mean
+use_mel_loss: true   # Whether to use Mel-spectrogram loss.
+mel_loss_params:     # Mel-spectrogram loss parameters.
+    fs: 16000
+    fft_size: 1024
+    hop_size: 256
+    win_length: null
+    window: "hann"
+    num_mels: 80
+    fmin: 0
+    fmax: 8000
+    log_base: null   # Log base. If set to null, use natural logarithm.
+generator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+discriminator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+use_feat_match_loss: true
+feat_match_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+    average_by_layers: false         # Whether to average loss by #layers in each discriminator.
+    include_final_outputs: true      # Whether to include final outputs in feat match loss calculation.
+
+###########################################################
+#               ADVERSARIAL LOSS SETTING                  #
+###########################################################
+lambda_aux: 45.0       # Loss balancing coefficient for STFT loss.
+lambda_adv: 1.0        # Loss balancing coefficient for adversarial loss.
+lambda_feat_match: 2.0 # Loss balancing coefficient for feat match loss..
+
+###########################################################
+#                  DATA LOADER SETTING                    #
+###########################################################
+batch_size: 16              # Batch size.
+batch_max_steps: 10240      # Length of each audio in batch. Make sure dividable by hop_size.
+pin_memory: true            # Whether to pin memory in Pytorch DataLoader.
+num_workers: 0              # Number of workers in Pytorch DataLoader.
+remove_short_samples: false # Whether to remove samples the length of which are less than batch_max_steps.
+allow_cache: true           # Whether to allow cache in dataset. If true, it requires cpu memory.
+
+###########################################################
+#             OPTIMIZER & SCHEDULER SETTING               #
+###########################################################
+generator_optimizer_type: Adam
+generator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+generator_scheduler_type: MultiStepLR
+generator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+generator_grad_norm: -1
+discriminator_optimizer_type: Adam
+discriminator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+discriminator_scheduler_type: MultiStepLR
+discriminator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+discriminator_grad_norm: -1
+
+###########################################################
+#                    INTERVAL SETTING                     #
+###########################################################
+generator_train_start_steps: 1     # Number of steps to start to train discriminator.
+discriminator_train_start_steps: 0 # Number of steps to start to train discriminator.
+train_max_steps: 400000           # Number of training steps.
+save_interval_steps: 50000         # Interval steps to save checkpoint.
+eval_interval_steps: 10000          # Interval steps to evaluate the network.
+log_interval_steps: 100            # Interval steps to record the training log.
+
+###########################################################
+#                     OTHER SETTING                       #
+###########################################################
+num_save_intermediate_results: 4  # Number of results to be saved as intermediate results.
diff --git a/egs/opencpop/token_voc1/conf/hifigan_token_duration_24k.v1.yaml b/egs/opencpop/token_voc1/conf/hifigan_token_duration_24k.v1.yaml
new file mode 100644
index 00000000..af1538a6
--- /dev/null
+++ b/egs/opencpop/token_voc1/conf/hifigan_token_duration_24k.v1.yaml
@@ -0,0 +1,186 @@
+# This configuration is based on HiFiGAN V1, derived
+# from official repository (https://github.com/jik876/hifi-gan).
+
+###########################################################
+#                FEATURE EXTRACTION SETTING               #
+###########################################################
+sampling_rate: 24000     # Sampling rate.
+fft_size: null           # FFT size.
+hop_size: 320            # Hop size.
+win_length: null         # Window length.
+                         # If set to null, it will be the same as fft_size.
+window: null             # Window function.
+num_mels: 1              # Number of mel basis.
+fmin: null               # Minimum freq in mel basis calculation.
+fmax: null               # Maximum frequency in mel basis calculation.
+global_gain_scale: 1.0   # Will be multiplied to all of waveform.
+trim_silence: false      # Whether to trim the start and end of silence.
+trim_threshold_in_db: 20 # Need to tune carefully if the recording is not good.
+trim_frame_size: 1024    # Frame size in trimming.
+trim_hop_size: 256       # Hop size in trimming.
+format: "hdf5"           # Feature file format. "npy" or "hdf5" is supported.
+
+###########################################################
+#         GENERATOR NETWORK ARCHITECTURE SETTING          #
+###########################################################
+generator_type: DiscreteSymbolF0Generator
+generator_params:
+    in_channels: 768                      # Number of input channels.
+    out_channels: 1                       # Number of output channels.
+    channels: 512                         # Number of initial channels.
+    linear_channel: 256                   # Number of Linear channels(f0).
+    num_embs: 1025
+    kernel_size: 7                        # Kernel size of initial and final conv layers.
+    upsample_scales: [10, 8, 2, 2]        # Upsampling scales.
+    upsample_kernel_sizes: [20, 16, 4, 4] # Kernel size for upsampling layers.
+    resblock_kernel_sizes: [3, 7, 11]     # Kernel size for residual blocks.
+    resblock_dilations:                   # Dilations for residual blocks.
+        - [1, 3, 5]
+        - [1, 3, 5]
+        - [1, 3, 5]
+    use_additional_convs: true            # Whether to use additional conv layer in residual blocks.
+    bias: true                            # Whether to use bias parameter in conv.
+    nonlinear_activation: "LeakyReLU"     # Nonlinear activation type.
+    nonlinear_activation_params:          # Nonlinear activation paramters.
+        negative_slope: 0.1
+    use_weight_norm: true                 # Whether to apply weight normalization.
+    # duration_layers: 2                    # Duration predictor layers
+    # duration_chans: 384                   # Duration predictor conv channels
+    # duration_kernel_size: 3               # Duration predictor kernel size
+    # duration_offset: 1.0                  # Duration predictor offset 
+    # duration_dropout_rate: 0.5            # Duration predictor dropout
+    num_spk_embs: 0                       # Do not consider speaker embedding for single spk
+    
+    use_weight_sum: false                  # Whether to apply weight sum
+
+###########################################################
+#       DISCRIMINATOR NETWORK ARCHITECTURE SETTING        #
+###########################################################
+discriminator_type: HiFiGANMultiScaleMultiPeriodDiscriminator
+discriminator_params:
+    scales: 3                              # Number of multi-scale discriminator.
+    scale_downsample_pooling: "AvgPool1d"  # Pooling operation for scale discriminator.
+    scale_downsample_pooling_params:
+        kernel_size: 4                     # Pooling kernel size.
+        stride: 2                          # Pooling stride.
+        padding: 2                         # Padding size.
+    scale_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [15, 41, 5, 3]       # List of kernel sizes.
+        channels: 128                      # Initial number of channels.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        max_groups: 16                     # Maximum number of groups in downsampling conv layers.
+        bias: true
+        downsample_scales: [4, 4, 4, 4, 1] # Downsampling scales.
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:
+            negative_slope: 0.1
+    follow_official_norm: true             # Whether to follow the official norm setting.
+    periods: [2, 3, 5, 7, 11]              # List of period for multi-period discriminator.
+    period_discriminator_params:
+        in_channels: 1                     # Number of input channels.
+        out_channels: 1                    # Number of output channels.
+        kernel_sizes: [5, 3]               # List of kernel sizes.
+        channels: 32                       # Initial number of channels.
+        downsample_scales: [3, 3, 3, 3, 1] # Downsampling scales.
+        max_downsample_channels: 1024      # Maximum number of channels in downsampling conv layers.
+        bias: true                         # Whether to use bias parameter in conv layer."
+        nonlinear_activation: "LeakyReLU"  # Nonlinear activation.
+        nonlinear_activation_params:       # Nonlinear activation paramters.
+            negative_slope: 0.1
+        use_weight_norm: true              # Whether to apply weight normalization.
+        use_spectral_norm: false           # Whether to apply spectral normalization.
+
+###########################################################
+#                   STFT LOSS SETTING                     #
+###########################################################
+use_stft_loss: false # Whether to use multi-resolution STFT loss.
+use_duration_loss: true # Whether to use duration prediction loss (need for duration predictor)
+duration_loss_params:
+    offset: 1.0
+    reduction: mean
+use_mel_loss: true   # Whether to use Mel-spectrogram loss.
+mel_loss_params:     # Mel-spectrogram loss parameters.
+    fs: 24000
+    fft_size: 2048
+    hop_size: 300
+    win_length: 1200
+    window: "hann"
+    num_mels: 80
+    fmin: 0
+    fmax: 8000
+    log_base: null   # Log base. If set to null, use natural logarithm.
+generator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+discriminator_adv_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+use_feat_match_loss: true
+feat_match_loss_params:
+    average_by_discriminators: false # Whether to average loss by #discriminators.
+    average_by_layers: false         # Whether to average loss by #layers in each discriminator.
+    include_final_outputs: true      # Whether to include final outputs in feat match loss calculation.
+
+###########################################################
+#               ADVERSARIAL LOSS SETTING                  #
+###########################################################
+lambda_aux: 45.0       # Loss balancing coefficient for STFT loss.
+lambda_adv: 1.0        # Loss balancing coefficient for adversarial loss.
+lambda_feat_match: 2.0 # Loss balancing coefficient for feat match loss..
+
+###########################################################
+#                  DATA LOADER SETTING                    #
+###########################################################
+batch_size: 16              # Batch size.
+batch_max_steps: 10240      # Length of each audio in batch. Make sure dividable by hop_size.
+pin_memory: true            # Whether to pin memory in Pytorch DataLoader.
+num_workers: 0              # Number of workers in Pytorch DataLoader.
+remove_short_samples: false # Whether to remove samples the length of which are less than batch_max_steps.
+allow_cache: true           # Whether to allow cache in dataset. If true, it requires cpu memory.
+
+###########################################################
+#             OPTIMIZER & SCHEDULER SETTING               #
+###########################################################
+generator_optimizer_type: Adam
+generator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+generator_scheduler_type: MultiStepLR
+generator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+generator_grad_norm: -1
+discriminator_optimizer_type: Adam
+discriminator_optimizer_params:
+    lr: 2.0e-4
+    betas: [0.5, 0.9]
+    weight_decay: 0.0
+discriminator_scheduler_type: MultiStepLR
+discriminator_scheduler_params:
+    gamma: 0.5
+    milestones:
+        - 200000
+        - 400000
+        - 600000
+        - 800000
+discriminator_grad_norm: -1
+
+###########################################################
+#                    INTERVAL SETTING                     #
+###########################################################
+generator_train_start_steps: 1     # Number of steps to start to train discriminator.
+discriminator_train_start_steps: 0 # Number of steps to start to train discriminator.
+train_max_steps: 250000           # Number of training steps.
+save_interval_steps: 50000         # Interval steps to save checkpoint.
+eval_interval_steps: 10000          # Interval steps to evaluate the network.
+log_interval_steps: 100            # Interval steps to record the training log.
+
+###########################################################
+#                     OTHER SETTING                       #
+###########################################################
+num_save_intermediate_results: 4  # Number of results to be saved as intermediate results.
diff --git a/egs/opencpop/token_voc1/conf/slurm.conf b/egs/opencpop/token_voc1/conf/slurm.conf
new file mode 100644
index 00000000..cefd21f0
--- /dev/null
+++ b/egs/opencpop/token_voc1/conf/slurm.conf
@@ -0,0 +1,12 @@
+# Default configuration
+command sbatch --export=PATH  --ntasks-per-node=1
+option time=* --time $0
+option mem=* --mem-per-cpu $0
+option mem=0          # Do not add anything to qsub_opts
+option num_threads=* --cpus-per-task $0 --ntasks-per-node=1
+option num_threads=1 --cpus-per-task 1  --ntasks-per-node=1 # Do not add anything to qsub_opts
+default gpu=0
+option gpu=0 -p cpu
+option gpu=* -p gpu --gres=gpu:$0
+# note: the --max-jobs-run option is supported as a special case
+# by slurm.pl and you don't have to handle it in the config file.
diff --git a/egs/opencpop/token_voc1/local/data_prep.py b/egs/opencpop/token_voc1/local/data_prep.py
new file mode 100644
index 00000000..cf3ce9fb
--- /dev/null
+++ b/egs/opencpop/token_voc1/local/data_prep.py
@@ -0,0 +1,62 @@
+import argparse
+import os
+
+import librosa
+import numpy as np
+
+
+def process_utterance(wavscp, utt2spk, audio_dir, wav_dumpdir, segment, tgt_sr=24000):
+    uid, _, phns, notes, syb_dur, phn_dur, keep = segment.strip().split("|")
+
+    utt2spk.write("{} {}\n".format(uid, "opencpop"))
+
+    # apply bit convert, there is a known issue in direct convert in format wavscp
+    cmd = (
+        f"sox {os.path.join(audio_dir, uid)}.wav -c 1 -t wavpcm -b 16 -r"
+        f" {tgt_sr} {os.path.join(wav_dumpdir, uid)}.wav"
+    )
+    os.system(cmd)
+
+    wavscp.write("{} {}.wav\n".format(uid, os.path.join(wav_dumpdir, uid)))
+
+
+def process_subset(args, set_name):
+    if not os.path.exists(args.tgt_dir):
+        os.makedirs(args.tgt_dir)
+    if not os.path.exists(os.path.join(args.tgt_dir, set_name)):
+        os.makedirs(os.path.join(args.tgt_dir, set_name))
+    wavscp = open(
+        os.path.join(args.tgt_dir, set_name, "wav.scp"), "w", encoding="utf-8"
+    )
+    utt2spk = open(
+        os.path.join(args.tgt_dir, set_name, "utt2spk"), "w", encoding="utf-8"
+    )
+    with open(
+        os.path.join(args.src_data, "segments", set_name + ".txt"),
+        "r",
+        encoding="utf-8",
+    ) as f:
+        segments = f.read().strip().split("\n")
+        for segment in segments:
+            process_utterance(
+                wavscp,
+                utt2spk,
+                os.path.join(args.src_data, "segments", "wavs"),
+                args.wav_dumpdir,
+                segment,
+                tgt_sr=args.sr,
+            )
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Prepare Data for Opencpop Database")
+    parser.add_argument("src_data", type=str, help="source data directory")
+    parser.add_argument("--tgt_dir", type=str, default="data")
+    parser.add_argument(
+        "--wav_dumpdir", type=str, help="wav dump directory (rebit)", default="wav_dump"
+    )
+    parser.add_argument("--sr", type=int, help="sampling rate (Hz)")
+    args = parser.parse_args()
+
+    for name in ["train", "test"]:
+        process_subset(args, name)
diff --git a/egs/opencpop/token_voc1/local/preprocess_token.py b/egs/opencpop/token_voc1/local/preprocess_token.py
new file mode 100755
index 00000000..75a2f65c
--- /dev/null
+++ b/egs/opencpop/token_voc1/local/preprocess_token.py
@@ -0,0 +1,413 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2019 Tomoki Hayashi
+#  MIT License (https://opensource.org/licenses/MIT)
+
+"""Perform preprocessing and raw feature extraction."""
+
+import argparse
+import logging
+import os
+
+import librosa
+import numpy as np
+import resampy
+import soundfile as sf
+import yaml
+import torch
+from tqdm import tqdm
+from scipy.interpolate import interp1d
+
+from parallel_wavegan.datasets import AudioDataset, AudioSCPDataset
+from parallel_wavegan.utils import write_hdf5
+
+
+def _convert_to_continuous_f0(f0: np.array) -> np.array:
+    if (f0 == 0).all():
+        logging.warning("All frames seems to be unvoiced.")
+        return f0
+
+    # padding start and end of f0 sequence
+    start_f0 = f0[f0 != 0][0]
+    end_f0 = f0[f0 != 0][-1]
+    start_idx = np.where(f0 == start_f0)[0][0]
+    end_idx = np.where(f0 == end_f0)[0][-1]
+    f0[:start_idx] = start_f0
+    f0[end_idx:] = end_f0
+
+    # get non-zero frame index
+    nonzero_idxs = np.where(f0 != 0)[0]
+
+    # perform linear interpolation
+    interp_fn = interp1d(nonzero_idxs, f0[nonzero_idxs])
+    f0 = interp_fn(np.arange(0, f0.shape[0]))
+
+    return f0
+
+def f0_dio(
+    audio,
+    sampling_rate,
+    hop_size=160,
+    pitch_min=80,
+    pitch_max=10000,
+    use_log_f0=True,
+    use_continuous_f0=True,
+):
+    """Compute F0 with pyworld.dio
+
+    Args:
+        audio (ndarray): Audio signal (T,).
+        sampling_rate (int): Sampling rate.
+        hop_size (int): Hop size.
+        pitch_min (int): Minimum pitch in pitch extraction.
+        pitch_max (int): Maximum pitch in pitch extraction.
+
+    Returns:
+        ndarray: f0 feature (#frames, ).
+
+    Note:
+        Unvoiced frame has value = 0.
+
+    """
+    if torch.is_tensor(audio):
+        x = audio.cpu().numpy().astype(np.double)
+    else:
+        x = audio.astype(np.double)
+    frame_period = 1000 * hop_size / sampling_rate
+    import pyworld
+    f0, timeaxis = pyworld.dio(
+        x,
+        sampling_rate,
+        f0_floor=pitch_min,
+        f0_ceil=pitch_max,
+        frame_period=frame_period,
+    )
+    f0 = pyworld.stonemask(x, f0, timeaxis, sampling_rate)
+    if use_continuous_f0:
+        f0 = _convert_to_continuous_f0(f0)
+    if use_log_f0:
+        nonzero_idxs = np.where(f0 != 0)[0]
+        f0[nonzero_idxs] = np.log(f0[nonzero_idxs])
+    return f0
+
+
+def main():
+    """Run preprocessing process."""
+    parser = argparse.ArgumentParser(
+        description=(
+            "Preprocess audio and then extract features (See detail in"
+            " parallel_wavegan/bin/preprocess.py)."
+        )
+    )
+    parser.add_argument(
+        "--wav-scp",
+        "--scp",
+        default=None,
+        type=str,
+        help="kaldi-style wav.scp file. you need to specify either scp or rootdir.",
+    )
+    parser.add_argument(
+        "--segments",
+        default=None,
+        type=str,
+        help=(
+            "kaldi-style segments file. if use, you must to specify both scp and"
+            " segments."
+        ),
+    )
+    parser.add_argument(
+        "--text",
+        default=None,
+        type=str,
+        help="kaldi-style text format hubert embedding index.",
+    )
+    parser.add_argument(
+        "--utt2spk",
+        default=None,
+        type=str,
+        help=(
+            "kaldi-style utt2spk file. If you want to add global conditionning with "
+            "speaker id, you need to specify this argument."
+        ),
+    )
+    parser.add_argument(
+        "--spk2idx",
+        default=None,
+        type=str,
+        help=(
+            "kaldi-style spk2idx file. If you want to add global conditionning with "
+            "speaker id, you need to specify this argument."
+        ),
+    )
+    parser.add_argument(
+        "--rootdir",
+        default=None,
+        type=str,
+        help=(
+            "directory including wav files. you need to specify either scp or rootdir."
+        ),
+    )
+    parser.add_argument(
+        "--dumpdir",
+        type=str,
+        required=True,
+        help="directory to dump feature files.",
+    )
+    parser.add_argument(
+        "--config",
+        type=str,
+        required=True,
+        help="yaml format configuration file.",
+    )
+    parser.add_argument(
+        "--verbose",
+        type=int,
+        default=1,
+        help="logging level. higher is more logging. (default=1)",
+    )
+    parser.add_argument(
+        "--use-f0",
+        default=False,
+        action="store_true",
+        help="whether to use f0 sequence.",
+    )
+    parser.add_argument(
+        "--use-embedding-feats",
+        default=False,
+        action="store_true",
+        help="whether to use pretrain model to get feature.",
+    )
+    parser.add_argument(
+        "--emb-layer",
+        type=int,
+        default=1,
+        help="logging level. higher is more logging. (default=1)",
+    )
+    parser.add_argument(
+        "--pretrained-model",
+        type=str,
+        default="facebook/hubert-base-ls960",
+        help="pretrained model for embedding feature",
+    )
+    args = parser.parse_args()
+
+    # set logger
+    if args.verbose > 1:
+        logging.basicConfig(
+            level=logging.DEBUG,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    elif args.verbose > 0:
+        logging.basicConfig(
+            level=logging.INFO,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    else:
+        logging.basicConfig(
+            level=logging.WARN,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+        logging.warning("Skip DEBUG/INFO messages")
+
+    # load config
+    with open(args.config) as f:
+        config = yaml.load(f, Loader=yaml.Loader)
+    config.update(vars(args))
+
+    # check arguments
+    if (args.wav_scp is not None and args.rootdir is not None) or (
+        args.wav_scp is None and args.rootdir is None
+    ):
+        raise ValueError("Please specify either --rootdir or --wav-scp.")
+
+    # get dataset
+    if args.rootdir is not None:
+        dataset = AudioDataset(
+            args.rootdir,
+            "*.wav",
+            audio_load_fn=sf.read,
+            return_utt_id=True,
+        )
+    else:
+        dataset = AudioSCPDataset(
+            args.wav_scp,
+            segments=args.segments,
+            return_utt_id=True,
+            return_sampling_rate=True,
+        )
+
+    # get token single layer / multi layer
+    if not os.path.isdir(args.text): # single layer token file
+        with open(args.text) as f:
+            lines = [line.strip() for line in f.readlines()]
+        text = {
+            line.split(maxsplit=1)[0]: line.split(maxsplit=1)[1].split() for line in lines
+        }
+    else:  # multi-stream: directory of token files
+        text = {}
+        for fname in os.listdir(args.text):
+            fpath = os.path.join(args.text, fname)
+            with open(fpath, 'r') as f:
+                lines = [line.strip() for line in f.readlines()]
+            for line in lines:
+                utt_name, tokens = line.split(maxsplit=1) # name, list
+                tokens = tokens.split() 
+                if text.get(utt_name) is None:
+                    text[utt_name] = []
+                text[utt_name].append(tokens)
+
+    # load spk2utt file
+    if args.utt2spk is not None:
+        with open(args.utt2spk) as f:
+            lines = [l.replace("\n", "") for l in f.readlines()]
+        utt2spk = {l.split()[0]: l.split()[1] for l in lines}
+        with open(args.spk2idx) as f:
+            lines = [l.replace("\n", "") for l in f.readlines()]
+        spk2idx = {l.split()[0]: int(l.split()[1]) for l in lines}
+
+    # check directly existence
+    if not os.path.exists(args.dumpdir):
+        os.makedirs(args.dumpdir, exist_ok=True)
+
+    # process each data
+    for utt_id, (audio, fs) in tqdm(dataset):
+        # check
+        assert len(audio.shape) == 1, f"{utt_id} seems to be multi-channel signal."
+        assert (
+            np.abs(audio).max() <= 1.0
+        ), f"{utt_id} seems to be different from 16 bit PCM."
+
+        # downsample
+        if fs != config["sampling_rate"]:
+            audio = resampy.resample(audio, fs, config["sampling_rate"], axis=0)
+
+        # trim silence
+        if config["trim_silence"]:
+            audio, _ = librosa.effects.trim(
+                audio,
+                top_db=config["trim_threshold_in_db"],
+                frame_length=config["trim_frame_size"],
+                hop_length=config["trim_hop_size"],
+            )
+        
+        # use feature embedding(for teacher-forcing)
+        if args.use_embedding_feats:
+            os.environ["http_proxy"] = "http://127.0.0.1:7890"
+            os.environ["https_proxy"] = "http://127.0.0.1:7890"
+            from transformers import AutoModel, Wav2Vec2FeatureExtractor
+            pretrained_model = args.pretrained_model
+            logging.info(f'model: {pretrained_model}')
+            model = AutoModel.from_pretrained(pretrained_model, cache_dir='/data3/tyx/pretrain_model')
+            processor = Wav2Vec2FeatureExtractor.from_pretrained(pretrained_model, cache_dir='/data3/tyx/pretrain_model') 
+            # model = AutoModel.from_pretrained(pretrained_model, cache_dir='/data3/tyx/pretrain_model')
+            # processor = Wav2Vec2CTCTokenizer.from_pretrained(pretrained_model, cache_dir='/data3/tyx/pretrain_model') 
+            inputs = processor(audio, sampling_rate=fs, return_tensors="pt")
+            outputs = model(**inputs, output_hidden_states=True)
+            features = outputs.hidden_states
+            mel = features[args.emb_layer].squeeze(0).detach().numpy()
+            # mel: (T, 1024) 
+        else:
+            # use hubert index instead of mel
+            mel = np.array(text[utt_id]).astype(np.int64)
+            if mel.ndim > 1: 
+                mel = mel.transpose(1, 0)
+            else:
+                mel = mel.reshape(-1, 1)
+            # mel: (T, 1)
+        # logging.info(f'mel({mel.shape})')
+        
+        if args.spk2idx is not None:
+            spk = utt2spk[utt_id]
+            if spk in spk2idx:
+                idx = spk2idx[spk]
+            else:
+                logging.warn(f"{spk} is unknown speaker.")
+                max_idx = max(spk2idx.values()) + 1
+                idx = max_idx
+
+            # concatenate with mel
+            idx = np.repeat(np.array(idx).reshape(1, 1), len(mel), axis=0)
+            mel = np.concatenate([mel, idx], axis=1)
+
+        # make sure the audio length and feature length are matched
+        logging.info(f"Mod: {len(audio) - len(mel) * config['hop_size']}")
+        if len(mel) * config["hop_size"] <= len(audio):
+            logging.warning(
+                "len(mel) * config['hop_size'] <= len(audio), may be errors"
+            )
+        # logging.info(f'audio: {len(audio)}')
+        # import math
+        # token_len = math.ceil(len(audio) / config["hop_size"])
+        # assert len(mel) == token_len
+        mel = mel[: len(audio) // config["hop_size"]]
+        audio = audio[: len(mel) * config["hop_size"]]
+        assert len(mel) * config["hop_size"] == len(audio)
+        
+        # use f0
+        if args.use_f0:
+            f0 = f0_dio(
+                audio,
+                sampling_rate=config["sampling_rate"],
+                hop_size=config["hop_size"],
+            ) # (#frames,) 
+            # logging.info(f'f0({f0.shape}): {f0}') 
+            if len(f0) > len(mel):
+                f0 = f0[: len(mel)]
+            else:
+                f0 = np.pad(f0, (0, len(mel) - len(f0)), mode="edge")
+
+        # apply global gain
+        if config["global_gain_scale"] > 0.0:
+            audio *= config["global_gain_scale"]
+        if np.abs(audio).max() >= 1.0:
+            logging.warn(
+                f"{utt_id} causes clipping. "
+                "it is better to re-consider global gain scale."
+            )
+            continue
+
+        # save
+        if config["format"] == "hdf5":
+            # logging.info(f'mel: {mel.shape} f0: {f0.shape}')
+            write_hdf5(
+                os.path.join(args.dumpdir, f"{utt_id}.h5"),
+                "wave",
+                audio.astype(np.float32),
+            )
+            write_hdf5(
+                os.path.join(args.dumpdir, f"{utt_id}.h5"),
+                "feats",
+                mel.astype(np.float32),
+            )
+            if args.use_f0:
+                write_hdf5(
+                    os.path.join(args.dumpdir, f"{utt_id}.h5"),
+                    "f0",
+                    f0.astype(np.float32),
+                )
+                
+        elif config["format"] == "npy":
+            np.save(
+                os.path.join(args.dumpdir, f"{utt_id}-wave.npy"),
+                audio.astype(np.float32),
+                allow_pickle=False,
+            )
+            np.save(
+                os.path.join(args.dumpdir, f"{utt_id}-feats.npy"),
+                mel.astype(np.float32),
+                allow_pickle=False,
+            )
+            if args.use_f0:
+                np.save(
+                    os.path.join(args.dumpdir, f"{utt_id}-f0.npy"),
+                    f0.astype(np.float32),
+                    allow_pickle=False,
+                )
+                
+        else:
+            raise ValueError("support only hdf5 or npy format.")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/egs/opencpop/token_voc1/path.sh b/egs/opencpop/token_voc1/path.sh
new file mode 100644
index 00000000..38c57bf2
--- /dev/null
+++ b/egs/opencpop/token_voc1/path.sh
@@ -0,0 +1,33 @@
+# cuda related
+# export CUDA_HOME=/usr/local/cuda-11.3
+# export LD_LIBRARY_PATH="${CUDA_HOME}/lib64:${LD_LIBRARY_PATH}"
+
+# path related
+# export PRJ_ROOT="${PWD}/../../.."
+# if [ -e "${PRJ_ROOT}/tools/venv/bin/activate" ]; then
+#     # shellcheck disable=SC1090
+#     . "${PRJ_ROOT}/tools/venv/bin/activate"
+# fi
+
+# python related
+export OMP_NUM_THREADS=1
+export PYTHONIOENCODING=UTF-8
+export MPL_BACKEND=Agg
+
+# check installation
+if ! command -v parallel-wavegan-train > /dev/null; then
+    echo "Error: It seems setup is not finished." >&2
+    echo "Error: Please setup your environment by following README.md" >&2
+    return 1
+fi
+if ! command -v jq > /dev/null; then
+    echo "Error: It seems jq is not installed." >&2
+    echo "Error: Please install via \`sudo apt-get install jq\`." >&2
+    echo "Error: If you do not have sudo, please download from https://stedolan.github.io/jq/download/." >&2
+    return 1
+fi
+if ! command -v yq > /dev/null; then
+    echo "Error: It seems yq is not installed." >&2
+    echo "Error: Please install via \`pip install yq\`." >&2
+    return 1
+fi
diff --git a/egs/opencpop/token_voc1/run.sh b/egs/opencpop/token_voc1/run.sh
new file mode 100755
index 00000000..7c42c992
--- /dev/null
+++ b/egs/opencpop/token_voc1/run.sh
@@ -0,0 +1,255 @@
+#!/bin/bash
+
+# Copyright 2020 Tomoki Hayashi
+#  MIT License (https://opensource.org/licenses/MIT)
+
+. ./cmd.sh || exit 1;
+. ./path.sh || exit 1;
+
+# basic settings
+stage=-1       # stage to start
+stop_stage=100 # stage to stop
+verbose=1      # verbosity level (lower is less info)
+n_gpus=1       # number of gpus in training
+n_jobs=8      # number of parallel jobs in feature extraction
+
+# NOTE(kan-bayashi): renamed to conf to avoid conflict in parse_options.sh
+conf=conf/hifigan_token_16k_nodp_f0.v1.yaml
+
+# directory path setting
+db_root=/data3/tyx/dataset/opencpop # direcotry including wavfiles (MODIFY BY YOURSELF)
+                          # each wav filename in the directory should be unique
+                          # e.g.
+                          # /path/to/database
+                          # ├── utt_1.wav
+                          # ├── utt_2.wav
+                          # │   ...
+                          # └── utt_N.wav
+dumpdir=dump           # directory to dump features
+
+# training related setting
+tag=""     # tag for directory to save model
+resume=""  # checkpoint path to resume training
+           # (e.g. <path>/<to>/checkpoint-10000steps.pkl)
+
+# decoding related setting
+checkpoint="" # checkpoint path to be used for decoding
+              # if not provided, the latest one will be used
+              # (e.g. <path>/<to>/checkpoint-400000steps.pkl)
+
+train_set="train"       # name of training data directory
+dev_set="dev"           # name of development data direcotry
+eval_set="test"         # name of evaluation data direcotry
+
+token_text=""
+use_f0=true                    # whether to add f0 
+use_embedding_feats=false      # whether to use pretrain feature as input
+pretrained_model="facebook/hubert-base-ls960"      # pre-trained model (confirm it on Huggingface)
+emb_layer=6
+fs=16000
+subexp="exp"
+
+# shellcheck disable=SC1091
+. utils/parse_options.sh || exit 1;
+
+set -euo pipefail
+
+if [ "${stage}" -le 0 ] && [ "${stop_stage}" -ge 0 ]; then
+    echo "Stage 0: Data preparation"
+    if [ ! -e "${db_root}" ]; then
+        echo "ERROR: Opencpop data does not exist."
+    	echo "ERROR: Please download https://wenet.org.cn/opencpop/download/ and locate it at ${download_dir}"
+    	exit 1
+    fi
+    echo "Please make sure fs=${fs} is right sample rate for model."
+    mkdir -p wav_dump
+    python local/data_prep.py ${db_root} \
+        --wav_dumpdir wav_dump \
+        --sr ${fs} \
+
+    sort -o data/train/wav.scp data/train/wav.scp
+
+    dev_num=50
+    train_num=$(( $(wc -l < data/train/wav.scp) - dev_num ))
+
+    mkdir -p data/${dev_set}
+    head -n $train_num data/${train_set}/wav.scp > data/${train_set}/wav.scp.tmp
+    tail -n $dev_num data/${train_set}/wav.scp > data/${dev_set}/wav.scp.tmp
+    mv data/${dev_set}/wav.scp.tmp data/${dev_set}/wav.scp
+    mv data/${train_set}/wav.scp.tmp data/${train_set}/wav.scp
+    
+    head -n $train_num data/${train_set}/utt2spk > data/${train_set}/utt2spk.tmp
+    tail -n $dev_num data/${train_set}/utt2spk > data/${dev_set}/utt2spk.tmp
+    mv data/${dev_set}/utt2spk.tmp data/${dev_set}/utt2spk
+    mv data/${train_set}/utt2spk.tmp data/${train_set}/utt2spk
+fi
+
+if [ "${stage}" -le 1 ] && [ "${stop_stage}" -ge 1 ]; then
+    echo "Stage 1: Feature extraction"
+    if [ ! -e "${token_text}" ]; then
+        echo "Valid --token_text is not provided. Please prepare it by yourself."        
+        echo "--token_text have 2 kinds of input: path of token_text file / path of token files directory."
+        cat << EOF
+---------------------------------
+token_text file: like kaldi-style text as follows:
+utt_id_1 0 0 0 0 1 1 1 1 2 2 2 2
+utt_id_2 0 0 0 0 0 0 3 3 3 3 3 3 5 5 5 5
+...
+----------------------------------
+token files directory: token_text files described above
+It will run in multi-stream way, training set should update in conf/hifigan_token_16k_nodp_f0.v1.yaml.
+token files directory format as follows:
+token_dir/
+    - token_file(layer1)
+    - token_file(layer2)
+    ....
+EOF
+        exit 1
+    fi
+    # extract raw features
+    pids=()
+    for name in "${train_set}" "${dev_set}" "${eval_set}"; do
+    (
+        [ ! -e "${dumpdir}/${name}/raw" ] && mkdir -p "${dumpdir}/${name}/raw"
+        echo "Feature extraction start. See the progress via ${dumpdir}/${name}/raw/preprocessing.*.log."
+        utils/make_subset_data.sh "data/${name}" "${n_jobs}" "${dumpdir}/${name}/raw"
+
+        _opts=
+        if [ ${use_f0} == true ]; then
+            _opts+="--use-f0 "
+        fi
+        if [ ${use_embedding_feats} == "true" ]; then
+            _opts+="--use-embedding-feats "
+            _opts+="--pretrained-model ${pretrained_model} "
+            _opts+="--emb-layer ${emb_layer} "
+        fi
+
+        ${train_cmd} JOB=1:${n_jobs} "${dumpdir}/${name}/raw/preprocessing.JOB.log" \
+            local/preprocess_token.py \
+                --config "${conf}" \
+                --scp "${dumpdir}/${name}/raw/wav.JOB.scp" \
+                --dumpdir "${dumpdir}/${name}/raw/dump.JOB" \
+                --text "${token_text}" \
+                --verbose "${verbose}" ${_opts}
+        echo "Successfully finished feature extraction of ${name} set."
+    ) &
+    pids+=($!)
+    done
+    i=0; for pid in "${pids[@]}"; do wait "${pid}" || ((++i)); done
+    [ "${i}" -gt 0 ] && echo "$0: ${i} background jobs are failed." && exit 1;
+    echo "Successfully finished feature extraction."
+fi
+
+if [ -z "${tag}" ]; then
+    expdir="${subexp}/${train_set}_opencpop_$(basename "${conf}" .yaml)"
+else
+    expdir="${subexp}/${train_set}_opencpop_${tag}"
+fi
+
+if [ "${stage}" -le 2 ] && [ "${stop_stage}" -ge 2 ]; then
+    echo "Stage 2: Network training"
+    [ ! -e "${expdir}" ] && mkdir -p "${expdir}"
+    if [ "${n_gpus}" -gt 1 ]; then
+        train="python -m parallel_wavegan.distributed.launch --nproc_per_node ${n_gpus} -c parallel-wavegan-train"
+    else
+        train="parallel-wavegan-train"
+    fi
+    _opts=
+    if [ ${use_f0} == true ]; then
+        _opts+="--use-f0 "
+    fi
+    # shellcheck disable=SC2012
+    resume="$(ls -dt "${expdir}"/*.pkl | head -1 || true)"
+    echo "Training start. See the progress via ${expdir}/train.log."
+    ${cuda_cmd} --gpu "${n_gpus}" "${expdir}/train.log" \
+        ${train} \
+            --config "${conf}" \
+            --train-dumpdir "${dumpdir}/${train_set}/raw" \
+            --dev-dumpdir "${dumpdir}/${dev_set}/raw" \
+            --outdir "${expdir}" \
+            --resume "${resume}" \
+            --verbose "${verbose}" ${_opts}
+    echo "Successfully finished training."
+fi
+
+if [ "${stage}" -le 3 ] && [ "${stop_stage}" -ge 3 ]; then
+    echo "Stage 3: Network decoding"
+    # shellcheck disable=SC2012
+    [ -z "${checkpoint}" ] && checkpoint="$(ls -dt "${expdir}"/*.pkl | head -1 || true)"
+    outdir="${expdir}/wav/$(basename "${checkpoint}" .pkl)"
+    pids=()
+    for name in "${dev_set}" "${eval_set}"; do
+    (
+        [ ! -e "${outdir}/${name}" ] && mkdir -p "${outdir}/${name}"
+        [ "${n_gpus}" -gt 1 ] && n_gpus=1
+        echo "Decoding start. See the progress via ${outdir}/${name}/decode.log."
+        _opts=
+        if [ ${use_f0} == true ]; then
+            _opts+="--use-f0 "
+        fi
+        ${cuda_cmd} --gpu "${n_gpus}" "${outdir}/${name}/decode.log" \
+            parallel-wavegan-decode \
+                --dumpdir "${dumpdir}/${name}/raw" \
+                --checkpoint "${checkpoint}" \
+                --outdir "${outdir}/${name}" \
+                --verbose "${verbose}" ${_opts}      
+        echo "Successfully finished decoding of ${name} set."
+    ) &
+    pids+=($!)
+    done
+    i=0; for pid in "${pids[@]}"; do wait "${pid}" || ((++i)); done
+    [ "${i}" -gt 0 ] && echo "$0: ${i} background jobs are failed." && exit 1;
+    echo "Successfully finished decoding."
+fi
+
+if [ "${stage}" -le 4 ] && [ "${stop_stage}" -ge 4 ]; then
+    echo "Stage 4: Scoring"
+    [ -z "${checkpoint}" ] && checkpoint="$(ls -dt "${expdir}"/*.pkl | head -1 || true)"
+    for dset in ${eval_set}; do
+        _data="data/${dset}"
+        _gt_wavscp="${_data}/wav.scp"
+        _dir="${expdir}/wav/$(basename "${checkpoint}" .pkl)"
+        _gen_wavdir="${_dir}/${dset}"
+
+        # Objective Evaluation - MCD
+        echo "Begin Scoring for MCD metrics on ${dset}, results are written under ${_dir}/MCD_res"
+
+        mkdir -p "${_dir}/MCD_res"
+        python utils/evaluate_mcd.py \
+            ${_gen_wavdir} \
+            ${_gt_wavscp} \
+            --outdir "${_dir}/MCD_res"
+
+        # Objective Evaluation - log-F0 RMSE
+        echo "Begin Scoring for F0 related metrics on ${dset}, results are written under ${_dir}/F0_res"
+
+        mkdir -p "${_dir}/F0_res"
+        python utils/evaluate_f0.py \
+            ${_gen_wavdir} \
+            ${_gt_wavscp} \
+            --outdir "${_dir}/F0_res"
+
+        # # Objective Evaluation - semitone ACC
+        # echo "Begin Scoring for SEMITONE related metrics on ${dset}, results are written under ${_dir}/SEMITONE_res"
+
+        # mkdir -p "${_dir}/SEMITONE_res"
+        # python utils/evaluate_semitone.py \
+        #     ${_gen_wavdir} \
+        #     ${_gt_wavscp} \
+        #     --outdir "${_dir}/SEMITONE_res"
+
+        #     # Objective Evaluation - VUV error
+        # echo "Begin Scoring for VUV related metrics on ${dset}, results are written under ${_dir}/VUV_res"
+
+        # mkdir -p "${_dir}/VUV_res"
+        # python utils/evaluate_vuv.py \
+        #     ${_gen_wavdir} \
+        #     ${_gt_wavscp} \
+        #     --outdir "${_dir}/VUV_res"
+
+    done
+else
+    echo "Skip the evaluation stages"
+fi
+
+echo "Finished."
diff --git a/egs/opencpop/token_voc1/utils b/egs/opencpop/token_voc1/utils
new file mode 120000
index 00000000..973afe67
--- /dev/null
+++ b/egs/opencpop/token_voc1/utils
@@ -0,0 +1 @@
+../../../utils
\ No newline at end of file
diff --git a/egs/opencpop/voc1/path.sh b/egs/opencpop/voc1/path.sh
index 393220c4..bf066c3f 100644
--- a/egs/opencpop/voc1/path.sh
+++ b/egs/opencpop/voc1/path.sh
@@ -1,13 +1,13 @@
-# cuda related
-export CUDA_HOME=/usr/local/cuda-11.3
-export LD_LIBRARY_PATH="${CUDA_HOME}/lib64:${LD_LIBRARY_PATH}"
+# # cuda related
+# export CUDA_HOME=/usr/local/cuda-11.3
+# export LD_LIBRARY_PATH="${CUDA_HOME}/lib64:${LD_LIBRARY_PATH}"
 
-# path related
-export PRJ_ROOT="${PWD}/../../../../Muskits/"
-if [ -e "${PRJ_ROOT}/tools/activate_python.sh" ]; then
-    # shellcheck disable=SC1090
-    . "${PRJ_ROOT}/tools/activate_python.sh"
-fi
+# # path related
+# export PRJ_ROOT="${PWD}/../../../../Muskits/"
+# if [ -e "${PRJ_ROOT}/tools/activate_python.sh" ]; then
+#     # shellcheck disable=SC1090
+#     . "${PRJ_ROOT}/tools/activate_python.sh"
+# fi
 
 # python related
 export OMP_NUM_THREADS=1
diff --git a/egs/opencpop/voc1/run.sh b/egs/opencpop/voc1/run.sh
index 920a51bb..c0caf52b 100755
--- a/egs/opencpop/voc1/run.sh
+++ b/egs/opencpop/voc1/run.sh
@@ -17,7 +17,7 @@ n_jobs=8      # number of parallel jobs in feature extraction
 conf=conf/uhifigan.v1.yaml
 
 # directory path setting
-download_dir= # set the directory to your database
+download_dir=/data3/tyx/dataset # set the directory to your database
 dumpdir=dump           # directory to dump features
 
 # training related setting
@@ -51,7 +51,7 @@ fi
 if [ "${stage}" -le 0 ] && [ "${stop_stage}" -ge 0 ]; then
     echo "Stage 0: Data preparation"
     mkdir -p wav_dump
-    python local/data_prep.py ${download_dir}/Opencpop \
+    python local/data_prep.py ${download_dir}/opencpop \
         --wav_dumpdir wav_dump \
         --sr 24000
 
@@ -81,7 +81,8 @@ if [ "${stage}" -le 1 ] && [ "${stop_stage}" -ge 1 ]; then
                 --config "${conf}" \
                 --scp "${dumpdir}/${name}/raw/wav.JOB.scp" \
                 --dumpdir "${dumpdir}/${name}/raw/dump.JOB" \
-                --verbose "${verbose}"
+                --verbose "${verbose}" \
+                --extract-f0
         echo "Successfully finished feature extraction of ${name} set."
     ) &
     pids+=($!)
@@ -97,7 +98,7 @@ if [ "${stage}" -le 1 ] && [ "${stop_stage}" -ge 1 ]; then
             --config "${conf}" \
             --rootdir "${dumpdir}/${train_set}/raw" \
             --dumpdir "${dumpdir}/${train_set}" \
-            --verbose "${verbose}"
+            --verbose "${verbose}" 
     echo "Successfully finished calculation of statistics."
 
     # normalize and dump them
diff --git a/parallel_wavegan/bin/decode.py b/parallel_wavegan/bin/decode.py
index 3c88bf08..9da5067a 100755
--- a/parallel_wavegan/bin/decode.py
+++ b/parallel_wavegan/bin/decode.py
@@ -23,6 +23,7 @@
     MelDataset,
     MelF0ExcitationDataset,
     MelSCPDataset,
+    MelF0Dataset,
 )
 from parallel_wavegan.utils import load_model, read_hdf5
 
@@ -96,6 +97,12 @@ def main():
         default=1,
         help="logging level. higher is more logging. (default=1)",
     )
+    parser.add_argument(
+        "--use-f0",
+        default=False,
+        action="store_true",
+        help="whether to use f0 sequence.",
+    )
     args = parser.parse_args()
 
     # set logger
@@ -164,6 +171,9 @@ def main():
             if config["format"] == "hdf5":
                 mel_query = "*.h5"
                 mel_load_fn = lambda x: read_hdf5(x, "feats")  # NOQA
+                if args.use_f0:
+                    f0_query = "*.h5"
+                    f0_load_fn = lambda x: read_hdf5(x, "f0")  # NOQA
                 if use_f0_and_excitation:
                     f0_query = "*.h5"
                     f0_load_fn = lambda x: read_hdf5(x, "f0")  # NOQA
@@ -172,6 +182,9 @@ def main():
             elif config["format"] == "npy":
                 mel_query = "*-feats.npy"
                 mel_load_fn = np.load
+                if args.use_f0:
+                    f0_query = "*-f0.npy"
+                    f0_load_fn = np.load
                 if use_f0_and_excitation:
                     f0_query = "*-f0.npy"
                     f0_load_fn = np.load
@@ -180,7 +193,16 @@ def main():
             else:
                 raise ValueError("Support only hdf5 or npy format.")
 
-            if not use_f0_and_excitation:
+            if args.use_f0:
+                dataset = MelF0Dataset(
+                    root_dir=args.dumpdir,
+                    mel_query=mel_query,
+                    f0_query=f0_query,
+                    mel_load_fn=mel_load_fn,
+                    f0_load_fn=f0_load_fn,
+                    return_utt_id=True,
+                )
+            elif not use_f0_and_excitation:
                 dataset = MelDataset(
                     args.dumpdir,
                     mel_query=mel_query,
@@ -213,7 +235,10 @@ def main():
         total_rtf = 0.0
         with torch.no_grad(), tqdm(dataset, desc="[decode]") as pbar:
             for idx, items in enumerate(pbar, 1):
-                if not use_f0_and_excitation:
+                if args.use_f0:
+                    utt_id, c, f0 = items
+                    excitation = None
+                elif not use_f0_and_excitation:
                     utt_id, c = items
                     f0, excitation = None, None
                 else:
diff --git a/parallel_wavegan/bin/train.py b/parallel_wavegan/bin/train.py
index 5fab384a..9e16dbb3 100755
--- a/parallel_wavegan/bin/train.py
+++ b/parallel_wavegan/bin/train.py
@@ -27,6 +27,7 @@
 from parallel_wavegan.datasets import (
     AudioDataset,
     AudioMelDataset,
+    AudioMelF0Dataset,
     AudioMelF0ExcitationDataset,
     AudioMelSCPDataset,
     AudioSCPDataset,
@@ -94,6 +95,9 @@ def __init__(
         self.use_duration_prediction = "Duration" in config.get(
             "generator_type", "ParallelWaveGANGenerator"
         )
+        # self.use_f0 = "F0" in config.get(
+        #     "generator_type", "ParallelWaveGANGenerator"
+        # )
 
     def run(self):
         """Run training."""
@@ -577,7 +581,7 @@ def _parse_batch(self, batch):
         """Parse batch and send to the device."""
         # parse batch
         if self.use_duration_prediction:
-            inputs, targets, durations = batch
+                inputs, targets, durations = batch
         else:
             inputs, targets = batch
 
@@ -652,6 +656,7 @@ def __init__(
         hop_size=256,
         aux_context_window=2,
         use_noise_input=False,
+        use_f0=False,
         use_f0_and_excitation=False,
         use_aux_input=True,
         use_duration=False,
@@ -682,6 +687,7 @@ def __init__(
         self.batch_max_steps = batch_max_steps
         self.aux_context_window = aux_context_window
         self.use_noise_input = use_noise_input
+        self.use_f0 = use_f0
         self.use_f0_and_excitation = use_f0_and_excitation
         self.use_aux_input = use_aux_input
         self.use_duration = use_duration
@@ -733,6 +739,8 @@ def __call__(self, batch):
                 self._adjust_length(*b) for b in batch if len(b[1]) > self.mel_threshold
             ]
             xs, cs = [b[0] for b in batch], [b[1] for b in batch]
+            if self.use_f0:
+                fs = [b[2] for b in batch]
             if self.use_f0_and_excitation:
                 fs, es = [b[2] for b in batch], [b[3] for b in batch]
 
@@ -764,6 +772,13 @@ def __call__(self, batch):
                 )  # (B, 1, T')
                 e_batch = torch.tensor(e_batch, dtype=torch.float)  # (B, 1, T', C')
                 e_batch = e_batch.reshape(e_batch.shape[0], 1, -1)  # (B, 1, T' * C')
+                
+            if self.use_f0:
+                f_batch = [f[start:end] for f, start, end in zip(fs, c_starts, c_ends)]
+                f_batch = np.array(f_batch)
+                f_batch = torch.tensor(f_batch, dtype=torch.float).unsqueeze(
+                    1
+                )  # (B, 1, T')
 
             # duration calculation and return with duration information
             if self.use_duration:
@@ -794,6 +809,8 @@ def __call__(self, batch):
                 input_items = (z_batch,) + input_items
             if self.use_f0_and_excitation:
                 input_items = input_items + (f_batch, e_batch)
+            if self.use_f0:
+                input_items = input_items + (f_batch,)
 
             return input_items, y_batch
         else:
@@ -892,6 +909,8 @@ def _adjust_length(self, x, c, f0=None, excitation=None):
 
         if f0 is not None and excitation is not None:
             return x, c, f0, excitation
+        elif f0 is not None and excitation is None:
+            return x, c, f0
         else:
             return x, c
 
@@ -1030,6 +1049,12 @@ def main():
         default=1,
         help="logging level. higher is more logging. (default=1)",
     )
+    parser.add_argument(
+        "--use-f0",
+        default=False,
+        action="store_true",
+        help="whether to use f0 sequence.",
+    )
     parser.add_argument(
         "--rank",
         "--local_rank",
@@ -1130,6 +1155,9 @@ def main():
                 f0_query, excitation_query = "*.h5", "*.h5"
                 f0_load_fn = lambda x: read_hdf5(x, "f0")  # NOQA
                 excitation_load_fn = lambda x: read_hdf5(x, "excitation")  # NOQA
+            if args.use_f0:
+                f0_query = "*.h5"
+                f0_load_fn = lambda x: read_hdf5(x, "f0")  # NOQA
             if use_local_condition:
                 local_query = "*.h5"
                 local_load_fn = lambda x: read_hdf5(x, "local")  # NOQA
@@ -1144,6 +1172,9 @@ def main():
                 f0_query, excitation_query = "*-f0.npy", "*-excitation.npy"
                 f0_load_fn = np.load
                 excitation_load_fn = np.load
+            if args.use_f0:
+                f0_query = "*-f0.npy"
+                f0_load_fn = np.load
             if use_local_condition:
                 local_query = "*-local.npy"
                 local_load_fn = np.load
@@ -1165,7 +1196,19 @@ def main():
 
     # define dataset for training data
     if args.train_dumpdir is not None:
-        if not use_f0_and_excitation:
+        if args.use_f0:
+            train_dataset = AudioMelF0Dataset(
+                root_dir=args.train_dumpdir,
+                audio_query=audio_query,
+                mel_query=mel_query,
+                f0_query=f0_query,
+                audio_load_fn=audio_load_fn,
+                mel_load_fn=mel_load_fn,
+                f0_load_fn=f0_load_fn,
+                mel_length_threshold=mel_length_threshold,
+                allow_cache=config.get("allow_cache", False),  # keep compatibility
+            )
+        elif not use_f0_and_excitation:
             if use_aux_input:
                 train_dataset = AudioMelDataset(
                     root_dir=args.train_dumpdir,
@@ -1211,6 +1254,10 @@ def main():
             raise NotImplementedError(
                 "SCP format is not supported for f0 and excitation."
             )
+        if args.use_f0:
+            raise NotImplementedError(
+                "SCP format is not supported for f0."
+            )
         if use_local_condition:
             raise NotImplementedError("Not supported.")
         if use_global_condition:
@@ -1233,7 +1280,19 @@ def main():
 
     # define dataset for validation
     if args.dev_dumpdir is not None:
-        if not use_f0_and_excitation:
+        if args.use_f0:
+                dev_dataset = AudioMelF0Dataset(
+                    root_dir=args.dev_dumpdir,
+                    audio_query=audio_query,
+                    mel_query=mel_query,
+                    f0_query=f0_query,
+                    audio_load_fn=audio_load_fn,
+                    mel_load_fn=mel_load_fn,
+                    f0_load_fn=f0_load_fn,
+                    mel_length_threshold=mel_length_threshold,
+                    allow_cache=config.get("allow_cache", False),  # keep compatibility
+                )
+        elif not use_f0_and_excitation:
             if use_aux_input:
                 dev_dataset = AudioMelDataset(
                     root_dir=args.dev_dumpdir,
@@ -1279,6 +1338,10 @@ def main():
             raise NotImplementedError(
                 "SCP format is not supported for f0 and excitation."
             )
+        if args.use_f0:
+            raise NotImplementedError(
+                "SCP format is not supported for f0."
+            )
         if use_local_condition:
             raise NotImplementedError("Not supported.")
         if use_global_condition:
@@ -1312,6 +1375,7 @@ def main():
         batch_max_steps=config["batch_max_steps"],
         hop_size=config.get("hop_size", None),
         aux_context_window=config["generator_params"].get("aux_context_window", 0),
+        use_f0=args.use_f0,
         use_f0_and_excitation=use_f0_and_excitation,
         use_noise_input=use_noise_input,
         use_aux_input=use_aux_input,
diff --git a/parallel_wavegan/datasets/audio_mel_dataset.py b/parallel_wavegan/datasets/audio_mel_dataset.py
index a375c51f..6adfffdb 100644
--- a/parallel_wavegan/datasets/audio_mel_dataset.py
+++ b/parallel_wavegan/datasets/audio_mel_dataset.py
@@ -191,6 +191,151 @@ def __len__(self):
         """
         return len(self.audio_files)
 
+class AudioMelF0Dataset(Dataset):
+    """PyTorch compatible audio and mel dataset."""
+
+    def __init__(
+        self,
+        root_dir,
+        audio_query="*.h5",
+        mel_query="*.h5",
+        f0_query="*.h5",
+        audio_load_fn=lambda x: read_hdf5(x, "wave"),
+        mel_load_fn=lambda x: read_hdf5(x, "feats"),
+        f0_load_fn=lambda x: read_hdf5(x, "f0"),
+        audio_length_threshold=None,
+        mel_length_threshold=None,
+        return_utt_id=False,
+        allow_cache=False,
+    ):
+        """Initialize dataset.
+
+        Args:
+            root_dir (str): Root directory including dumped files.
+            audio_query (str): Query to find audio files in root_dir.
+            mel_query (str): Query to find mel feature files in root_dir.
+            f0_query (str): Query to find f0 feature files in root_dir.
+            audio_load_fn (func): Function to load audio file.
+            mel_load_fn (func): Function to load mel feature file.
+            audio_load_fn (func): Function to load audio file.
+            mel_load_fn (func): Function to load mel feature file.
+            f0_load_fn (func): Function to load f0 feature file.
+            audio_length_threshold (int): Threshold to remove short audio files.
+            mel_length_threshold (int): Threshold to remove short feature files.
+            return_utt_id (bool): Whether to return the utterance id with arrays.
+            allow_cache (bool): Whether to allow cache of the loaded files.
+
+        """
+        # find all of audio and mel files
+        audio_files = sorted(find_files(root_dir, audio_query))
+        mel_files = sorted(find_files(root_dir, mel_query))
+        f0_files = sorted(find_files(root_dir, f0_query))
+
+        # filter by threshold
+        if audio_length_threshold is not None:
+            audio_lengths = [audio_load_fn(f).shape[0] for f in audio_files]
+            idxs = [
+                idx
+                for idx in range(len(audio_files))
+                if audio_lengths[idx] > audio_length_threshold
+            ]
+            if len(audio_files) != len(idxs):
+                logging.warning(
+                    "Some files are filtered by audio length threshold "
+                    f"({len(audio_files)} -> {len(idxs)})."
+                )
+            audio_files = [audio_files[idx] for idx in idxs]
+            mel_files = [mel_files[idx] for idx in idxs]
+            f0_files = [f0_files[idx] for idx in idxs]
+        if mel_length_threshold is not None:
+            mel_lengths = [mel_load_fn(f).shape[0] for f in mel_files]
+            idxs = [
+                idx
+                for idx in range(len(mel_files))
+                if mel_lengths[idx] > mel_length_threshold
+            ]
+            if len(mel_files) != len(idxs):
+                logging.warning(
+                    "Some files are filtered by mel length threshold "
+                    f"({len(mel_files)} -> {len(idxs)})."
+                )
+            audio_files = [audio_files[idx] for idx in idxs]
+            mel_files = [mel_files[idx] for idx in idxs]
+            f0_files = [f0_files[idx] for idx in idxs]
+
+        # assert the number of files
+        assert len(audio_files) != 0, f"Not found any audio files in ${root_dir}."
+        assert len(audio_files) == len(mel_files), (
+            f"Number of audio and mel files are different ({len(audio_files)} vs"
+            f" {len(mel_files)})."
+        )
+        assert len(audio_files) == len(f0_files), (
+            f"Number of audio and f0 files are different ({len(audio_files)} vs"
+            f" {len(f0_files)})."
+        )
+
+        self.audio_files = audio_files
+        self.audio_load_fn = audio_load_fn
+        self.mel_files = mel_files
+        self.mel_load_fn = mel_load_fn
+        self.f0_files = f0_files
+        self.f0_load_fn = f0_load_fn
+        if ".npy" in audio_query:
+            self.utt_ids = [
+                os.path.basename(f).replace("-wave.npy", "") for f in audio_files
+            ]
+        else:
+            self.utt_ids = [
+                os.path.splitext(os.path.basename(f))[0] for f in audio_files
+            ]
+
+        self.return_utt_id = return_utt_id
+        self.allow_cache = allow_cache
+        if allow_cache:
+            # NOTE(kan-bayashi): Manager is need to share memory in dataloader with num_workers > 0
+            self.manager = Manager()
+            self.caches = self.manager.list()
+            self.caches += [() for _ in range(len(audio_files))]
+
+    def __getitem__(self, idx):
+        """Get specified idx items.
+
+        Args:
+            idx (int): Index of the item.
+
+        Returns:
+            str: Utterance id (only in return_utt_id = True).
+            ndarray: Audio signal (T,).
+            ndarray: Feature (T', C).
+            ndarray: Feature (T', ).
+            ndarray: Feature (T', C').
+
+        """
+        if self.allow_cache and len(self.caches[idx]) != 0:
+            return self.caches[idx]
+
+        utt_id = self.utt_ids[idx]
+        audio = self.audio_load_fn(self.audio_files[idx])
+        mel = self.mel_load_fn(self.mel_files[idx])
+        f0 = self.f0_load_fn(self.f0_files[idx])
+
+        if self.return_utt_id:
+            items = utt_id, audio, mel, f0
+        else:
+            items = audio, mel, f0
+        if self.allow_cache:
+            self.caches[idx] = items
+        return items
+
+    def __len__(self):
+        """Return dataset length.
+
+        Returns:
+            int: The length of dataset.
+
+        """
+        return len(self.audio_files)
+
 
 class AudioMelF0ExcitationDataset(Dataset):
     """PyTorch compatible audio and mel dataset."""
@@ -755,3 +900,110 @@ def __len__(self):
 
         """
         return len(self.mel_files)
+
+
+class MelF0Dataset(Dataset):
+    """PyTorch compatible mel dataset."""
+
+    def __init__(
+        self,
+        root_dir,
+        mel_query="*-feats.npy",
+        f0_query="*-f0.npy",
+        mel_length_threshold=None,
+        mel_load_fn=np.load,
+        f0_load_fn=np.load,
+        return_utt_id=False,
+        allow_cache=False,
+    ):
+        """Initialize dataset.
+
+        Args:
+            root_dir (str): Root directory including dumped files.
+            mel_query (str): Query to find feature files in root_dir.
+            mel_load_fn (func): Function to load feature file.
+            mel_length_threshold (int): Threshold to remove short feature files.
+            return_utt_id (bool): Whether to return the utterance id with arrays.
+            allow_cache (bool): Whether to allow cache of the loaded files.
+
+        """
+        # find all of the mel files
+        mel_files = sorted(find_files(root_dir, mel_query))
+        f0_files = sorted(find_files(root_dir, f0_query))
+
+        # filter by threshold
+        if mel_length_threshold is not None:
+            mel_lengths = [mel_load_fn(f).shape[0] for f in mel_files]
+            idxs = [
+                idx
+                for idx in range(len(mel_files))
+                if mel_lengths[idx] > mel_length_threshold
+            ]
+            if len(mel_files) != len(idxs):
+                logging.warning(
+                    "Some files are filtered by mel length threshold "
+                    f"({len(mel_files)} -> {len(idxs)})."
+                )
+            mel_files = [mel_files[idx] for idx in idxs]
+            f0_files = [f0_files[idx] for idx in idxs]
+
+        # assert the number of files
+        assert len(mel_files) != 0, f"Not found any mel files in ${root_dir}."
+        assert len(f0_files) != 0, f"Not found any f0 files in ${root_dir}."
+
+        self.mel_files = mel_files
+        self.mel_load_fn = mel_load_fn
+        self.f0_files = f0_files
+        self.f0_load_fn = f0_load_fn
+
+        self.utt_ids = [os.path.splitext(os.path.basename(f))[0] for f in mel_files]
+        if ".npy" in mel_query:
+            self.utt_ids = [
+                os.path.basename(f).replace("-feats.npy", "") for f in mel_files
+            ]
+        else:
+            self.utt_ids = [os.path.splitext(os.path.basename(f))[0] for f in mel_files]
+        self.return_utt_id = return_utt_id
+        self.allow_cache = allow_cache
+        if allow_cache:
+            # NOTE(kan-bayashi): Manager is need to share memory in dataloader with num_workers > 0
+            self.manager = Manager()
+            self.caches = self.manager.list()
+            self.caches += [() for _ in range(len(mel_files))]
+
+    def __getitem__(self, idx):
+        """Get specified idx items.
+
+        Args:
+            idx (int): Index of the item.
+
+        Returns:
+            str: Utterance id (only in return_utt_id = True).
+            ndarray: Feature (T', C).
+
+        """
+        if self.allow_cache and len(self.caches[idx]) != 0:
+            return self.caches[idx]
+
+        utt_id = self.utt_ids[idx]
+        mel = self.mel_load_fn(self.mel_files[idx])
+        f0 = self.f0_load_fn(self.f0_files[idx])
+
+        if self.return_utt_id:
+            items = utt_id, mel, f0
+        else:
+            items = mel, f0
+
+        if self.allow_cache:
+            self.caches[idx] = items
+
+        return items
+
+    def __len__(self):
+        """Return dataset length.
+
+        Returns:
+            int: The length of dataset.
+
+        """
+        return len(self.mel_files)
diff --git a/parallel_wavegan/models/hifigan.py b/parallel_wavegan/models/hifigan.py
index 6dde0458..90dc31a2 100644
--- a/parallel_wavegan/models/hifigan.py
+++ b/parallel_wavegan/models/hifigan.py
@@ -886,6 +886,7 @@ def __init__(
         nonlinear_activation="LeakyReLU",
         nonlinear_activation_params={"negative_slope": 0.1},
         use_weight_norm=True,
+        use_embedding_feats=False,
     ):
         """Initialize HiFiGANGenerator module.
 
@@ -912,6 +913,7 @@ def __init__(
         """
         super().__init__()
         self.num_spk_embs = num_spk_embs
+        self.use_embedding_feats = use_embedding_feats
 
         # define id embedding
         self.emb = torch.nn.Embedding(
@@ -1015,9 +1017,11 @@ def forward(self, c):
             else:
                 g = g.unsqueeze(1).expand(-1, c.size(1), -1)
                 c = torch.cat([c, g], dim=-1)
+        # NOTE(Yuxun): update for using pretrain model layer output as input
         else:
-            assert c.size(1) == 1
-            c = self.emb(c.squeeze(1).long()).transpose(1, 2)  # (B, C, T)
+            if not self.use_embedding_feats:
+                assert c.size(1) == 1
+                c = self.emb(c.squeeze(1).long()).transpose(1, 2)  # (B, C, T)
 
         c = self.input_conv(c)
         for i in range(self.num_upsamples):
@@ -1085,8 +1089,10 @@ def inference(self, c, g=None, normalize_before=False):
         if g is not None:
             c = c[:, 0:1]
             c = torch.cat([c, c.new_zeros(*c.size()).fill_(g).to(c.device)], dim=1)
-        if self.num_spk_embs <= 0:
-            c = c[:, 0:1]
+        # NOTE(Yuxun): update for using pretrain model layer output as input  
+        # if self.num_spk_embs <= 0:
+        #     c = c[:, 0:1]
+        # c = self.forward(c.transpose(1, 0).unsqueeze(0))
         c = self.forward(c.transpose(1, 0).unsqueeze(0))
         return c.squeeze(0).transpose(1, 0)
 
@@ -1287,3 +1293,195 @@ def synthesis(self, c):
         c = self.output_conv(c)
 
         return c, ds_out
+
+
+class DiscreteSymbolF0Generator(DiscreteSymbolHiFiGANGenerator):
+    """Discrete Symbol HiFiGAN generator module with f0."""
+    
+    def __init__(
+        self,
+        in_channels=512,
+        out_channels=1,
+        channels=512,
+        linear_channel=256,
+        num_embs=100,
+        num_spk_embs=128,
+        spk_emb_dim=128,
+        concat_spk_emb=False,
+        kernel_size=7,
+        upsample_scales=(8, 8, 2, 2),
+        upsample_kernel_sizes=(16, 16, 4, 4),
+        resblock_kernel_sizes=(3, 7, 11),
+        resblock_dilations=[(1, 3, 5), (1, 3, 5), (1, 3, 5)],
+        use_additional_convs=True,
+        bias=True,
+        nonlinear_activation="LeakyReLU",
+        nonlinear_activation_params={"negative_slope": 0.1},
+        use_weight_norm=True,
+        # discret token
+        use_embedding_feats=False,
+        use_weight_sum=False,
+        layer_num=12,
+        use_fix_weight=False,
+        use_f0=True,
+    ):
+        """Initialize HiFiGANGenerator module.
+
+        Args:
+            in_channels (int): Number of input channels.
+            out_channels (int): Number of output channels.
+            channels (int): Number of hidden representation channels.
+            num_embs (int): Discrete symbol size
+            num_spk_embs (int): Speaker numbers for sPkeaer ID-based embedding
+            spk_emb_dim (int): Dimension of speaker embedding
+            concat_spk_emb (bool): whether to concat speaker embedding to the input
+            kernel_size (int): Kernel size of initial and final conv layer.
+            upsample_scales (list): List of upsampling scales.
+            upsample_kernel_sizes (list): List of kernel sizes for upsampling layers.
+            resblock_kernel_sizes (list): List of kernel sizes for residual blocks.
+            resblock_dilations (list): List of dilation list for residual blocks.
+            use_additional_convs (bool): Whether to use additional conv layers in residual blocks.
+            bias (bool): Whether to add bias parameter in convolution layers.
+            nonlinear_activation (str): Activation function module name.
+            nonlinear_activation_params (dict): Hyperparameters for activation function.
+            use_weight_norm (bool): Whether to use weight norm.
+                If set to true, it will be applied to all of the conv layers.
+
+        """
+        super().__init__(
+            in_channels=in_channels,
+            out_channels=out_channels,
+            channels=channels,
+            num_embs=num_embs,
+            num_spk_embs=num_spk_embs,
+            spk_emb_dim=spk_emb_dim,
+            concat_spk_emb=concat_spk_emb,
+            kernel_size=kernel_size,
+            upsample_scales=upsample_scales,
+            upsample_kernel_sizes=upsample_kernel_sizes,
+            resblock_kernel_sizes=resblock_kernel_sizes,
+            resblock_dilations=resblock_dilations,
+            use_additional_convs=use_additional_convs,
+            bias=bias,
+            nonlinear_activation=nonlinear_activation,
+            nonlinear_activation_params=nonlinear_activation_params,
+            use_weight_norm=use_weight_norm,
+        )
+        
+        if use_f0 is True:
+            self.use_f0 = use_f0
+            self.f0_embedding = torch.nn.Linear(
+                in_features=1,
+                out_features=linear_channel,
+            )
+        
+        self.use_weight_sum = use_weight_sum
+        if use_weight_sum is True:
+            self.layer_num = layer_num
+            self.weights = torch.nn.Parameter(torch.ones(self.layer_num))
+            self.use_fix_weight = use_fix_weight
+            
+            if use_fix_weight is True: # fix update
+                self.weights = torch.nn.Parameter(torch.ones(self.layer_num), requires_grad=False)
+            else: 
+                self.weights = torch.nn.Parameter(torch.ones(self.layer_num))
+                
+            self.emb = torch.nn.ModuleList([
+                torch.nn.Embedding(num_embeddings=num_embs, embedding_dim=in_channels) for _ in range(self.layer_num)
+            ])
+
+        self.input_conv = torch.nn.Conv1d(
+            in_channels + linear_channel if use_f0 is True else in_channels,
+            channels,
+            kernel_size,
+            1,
+            padding=(kernel_size - 1) // 2,
+        )
+
+        self.use_embedding_feats = use_embedding_feats
+        
+    
+    def forward(self, c, f0=None):
+        """Calculate forward propagation.
+
+        Args:
+            c (Tensor): Input tensor token: (B, 2, T). or (B, 1, T) or (B, L, T)
+                        or for embedding feature: (B, C, T)
+            f0 (Tensor): Input tensor (B, 1, T)
+        Returns:
+            Tensor: Output tensor (B, out_channels, T').
+        """
+        # logging.info(f'feats({c.shape}): {c}')       
+        # convert idx to embedding
+        if self.num_spk_embs > 0:
+            assert c.size(1) == 2
+            c_idx, g_idx = c.long().split(1, dim=1)
+            c = self.emb(c_idx.squeeze(1)).transpose(1, 2)  # (B, C, T)
+            g = self.spk_emb(g_idx[:, 0, 0])
+
+            # integrate global embedding
+            if not self.concat_spk_emb:
+                c = c + g.unsqueeze(2)
+            else:
+                g = g.unsqueeze(1).expand(-1, c.size(1), -1)
+                c = torch.cat([c, g], dim=-1)
+        else:
+            # NOTE(Yuxun): update for using pretrain model layer output as input
+            if not self.use_embedding_feats:
+                if self.use_weight_sum:
+                    assert c.size(1) == self.layer_num # (B, L, T)
+                    embedded = []
+                    for i, embedding_layer in enumerate(self.emb):
+                    # Apply the i-th embedding layer to the i-th layer of input
+                        embedded.append(embedding_layer(c[:, i].long()))
+                    c = torch.stack(embedded, dim=1).transpose(-1, 1)
+                    
+                    # weights: [L,]
+                    if self.use_fix_weight:
+                        norm_weights = self.weights
+                    else:
+                        norm_weights = F.softmax(self.weights, dim=-1) 
+                    # logging.info(f'norm_weights({norm_weights.shape}): {norm_weights}')
+                    # c: (B, C, T, L) * (L,) -> (B, C, T)
+                    c = torch.matmul(c, norm_weights)
+                else:
+                    assert c.size(1) == 1
+                    c = self.emb(c.squeeze(1).long()).transpose(1, 2)  # (B, C, T)
+         
+        # logging.info(f'f0({f0.shape}): {f0} ')        
+        if f0 is not None and self.use_f0:
+            f0 = self.f0_embedding(f0.transpose(1, 2)).transpose(1, 2)
+            c = torch.cat((c, f0), dim=1)
+        c = self.input_conv(c)
+        for i in range(self.num_upsamples):
+            c = self.upsamples[i](c)
+            cs = 0.0  # initialize
+            for j in range(self.num_blocks):
+                cs += self.blocks[i * self.num_blocks + j](c)
+            c = cs / self.num_blocks
+        c = self.output_conv(c)
+        return c
+    
+    
+    def inference(self, c, f0, g=None, normalize_before=False):
+        """Perform inference.
+
+        Args:
+            c (Union[Tensor, ndarray]): Input tensor (T, 2) or (T, 1) or (T, L).
+            f0 (Tensor): Input f0 (T,).
+        Returns:
+            Tensor: Output tensor (T ** prod(upsample_scales), out_channels).
+
+        """
+        assert not normalize_before, "No statistics are used."
+        if not isinstance(c, torch.Tensor):
+            c = torch.tensor(c, dtype=torch.long).to(next(self.parameters()).device)
+        if g is not None:
+            c = c[:, 0:1]
+            c = torch.cat([c, c.new_zeros(*c.size()).fill_(g).to(c.device)], dim=1)
+        if not self.use_embedding_feats:
+            if self.num_spk_embs <= 0 and not self.use_weight_sum:
+                c = c[:, 0:1]
+        # weight sum: c (T, L)
+        c = self.forward(c.transpose(1, 0).unsqueeze(0), f0.unsqueeze(0).unsqueeze(0))
+        return c.squeeze(0).transpose(1, 0)
\ No newline at end of file
diff --git a/utils/evaluate_f0.py b/utils/evaluate_f0.py
new file mode 100755
index 00000000..5464ed7d
--- /dev/null
+++ b/utils/evaluate_f0.py
@@ -0,0 +1,344 @@
+#!/usr/bin/env python3
+
+# Copyright 2021 Wen-Chin Huang and Tomoki Hayashi
+#  Apache 2.0  (http://www.apache.org/licenses/LICENSE-2.0)
+
+"""Evaluate log-F0 RMSE between generated and groundtruth audios based on World."""
+
+import argparse
+import fnmatch
+import logging
+import multiprocessing as mp
+import os
+from typing import Dict, List, Tuple
+
+import librosa
+import numpy as np
+import pysptk
+import pyworld as pw
+import soundfile as sf
+from fastdtw import fastdtw
+from scipy import spatial
+
+
+def find_files(
+    root_dir: str, query: List[str] = ["*.flac", "*.wav"], include_root_dir: bool = True
+) -> List[str]:
+    """Find files recursively.
+
+    Args:
+        root_dir (str): Root root_dir to find.
+        query (List[str]): Query to find.
+        include_root_dir (bool): If False, root_dir name is not included.
+
+    Returns:
+        List[str]: List of found filenames.
+
+    """
+    files = []
+    for root, dirnames, filenames in os.walk(root_dir, followlinks=True):
+        for q in query:
+            for filename in fnmatch.filter(filenames, q):
+                files.append(os.path.join(root, filename))
+    if not include_root_dir:
+        files = [file_.replace(root_dir + "/", "") for file_ in files]
+
+    return files
+
+
+def world_extract(
+    x: np.ndarray,
+    fs: int,
+    f0min: int = 40,
+    f0max: int = 800,
+    n_fft: int = 512,
+    n_shift: int = 256,
+    mcep_dim: int = 25,
+    mcep_alpha: float = 0.41,
+) -> np.ndarray:
+    """Extract World-based acoustic features.
+
+    Args:
+        x (ndarray): 1D waveform array.
+        fs (int): Minimum f0 value (default=40).
+        f0 (int): Maximum f0 value (default=800).
+        n_shift (int): Shift length in point (default=256).
+        n_fft (int): FFT length in point (default=512).
+        n_shift (int): Shift length in point (default=256).
+        mcep_dim (int): Dimension of mel-cepstrum (default=25).
+        mcep_alpha (float): All pass filter coefficient (default=0.41).
+
+    Returns:
+        ndarray: Mel-cepstrum with the size (N, n_fft).
+        ndarray: F0 sequence (N,).
+
+    """
+    # extract features
+    x = x.astype(np.float64)
+    f0, time_axis = pw.harvest(
+        x,
+        fs,
+        f0_floor=f0min,
+        f0_ceil=f0max,
+        frame_period=n_shift / fs * 1000,
+    )
+    sp = pw.cheaptrick(x, f0, time_axis, fs, fft_size=n_fft)
+    if mcep_dim is None or mcep_alpha is None:
+        mcep_dim, mcep_alpha = _get_best_mcep_params(fs)
+    mcep = pysptk.sp2mc(sp, mcep_dim, mcep_alpha)
+
+    return mcep, f0
+
+
+def _get_basename(path: str) -> str:
+    return os.path.splitext(os.path.split(path)[-1])[0]
+
+
+def _get_best_mcep_params(fs: int) -> Tuple[int, float]:
+    if fs == 16000:
+        return 23, 0.42
+    elif fs == 22050:
+        return 34, 0.45
+    elif fs == 24000:
+        return 34, 0.46
+    elif fs == 44100:
+        return 39, 0.53
+    elif fs == 48000:
+        return 39, 0.55
+    else:
+        raise ValueError(f"Not found the setting for {fs}.")
+
+
+def calculate(
+    file_list: List[str],
+    gt_file_list: List[str],
+    args: argparse.Namespace,
+    f0_rmse_dict: Dict[str, float],
+):
+    """Calculate log-F0 RMSE."""
+    for i, gen_path in enumerate(file_list):
+        corresponding_list = list(
+            filter(lambda gt_path: _get_basename(gt_path) in gen_path, gt_file_list)
+        )
+        assert len(corresponding_list) == 1
+        gt_path = corresponding_list[0]
+        gt_basename = _get_basename(gt_path)
+
+        # load wav file as int16
+        gen_x, gen_fs = sf.read(gen_path, dtype="int16")
+        gt_x, gt_fs = sf.read(gt_path, dtype="int16")
+
+        fs = gen_fs
+        if gen_fs != gt_fs:
+            gt_x = librosa.resample(gt_x.astype(np.float64), orig_sr=gt_fs, target_sr=gen_fs)
+
+        # extract ground truth and converted features
+        gen_mcep, gen_f0 = world_extract(
+            x=gen_x,
+            fs=fs,
+            f0min=args.f0min,
+            f0max=args.f0max,
+            n_fft=args.n_fft,
+            n_shift=args.n_shift,
+            mcep_dim=args.mcep_dim,
+            mcep_alpha=args.mcep_alpha,
+        )
+        gt_mcep, gt_f0 = world_extract(
+            x=gt_x,
+            fs=fs,
+            f0min=args.f0min,
+            f0max=args.f0max,
+            n_fft=args.n_fft,
+            n_shift=args.n_shift,
+            mcep_dim=args.mcep_dim,
+            mcep_alpha=args.mcep_alpha,
+        )
+
+        # DTW
+        _, path = fastdtw(gen_mcep, gt_mcep, dist=spatial.distance.euclidean)
+        twf = np.array(path).T
+        gen_f0_dtw = gen_f0[twf[0]]
+        gt_f0_dtw = gt_f0[twf[1]]
+
+        # Get voiced part
+        nonzero_idxs = np.where((gen_f0_dtw != 0) & (gt_f0_dtw != 0))[0]
+        gen_f0_dtw_voiced = np.log(gen_f0_dtw[nonzero_idxs])
+        gt_f0_dtw_voiced = np.log(gt_f0_dtw[nonzero_idxs])
+
+        # log F0 RMSE
+        log_f0_rmse = np.sqrt(np.mean((gen_f0_dtw_voiced - gt_f0_dtw_voiced) ** 2))
+        logging.info(f"{gt_basename} {log_f0_rmse:.4f}")
+        f0_rmse_dict[gt_basename] = log_f0_rmse
+
+
+def get_parser() -> argparse.Namespace:
+    """Get argument parser."""
+    parser = argparse.ArgumentParser(description="Evaluate Mel-cepstrum distortion.")
+    parser.add_argument(
+        "gen_wavdir_or_wavscp",
+        type=str,
+        help="Path of directory or wav.scp for generated waveforms.",
+    )
+    parser.add_argument(
+        "gt_wavdir_or_wavscp",
+        type=str,
+        help="Path of directory or wav.scp for ground truth waveforms.",
+    )
+    parser.add_argument(
+        "--outdir",
+        type=str,
+        help="Path of directory to write the results.",
+    )
+
+    # analysis related
+    parser.add_argument(
+        "--mcep_dim",
+        default=None,
+        type=int,
+        help=(
+            "Dimension of mel cepstrum coefficients. "
+            "If None, automatically set to the best dimension for the sampling."
+        ),
+    )
+    parser.add_argument(
+        "--mcep_alpha",
+        default=None,
+        type=float,
+        help=(
+            "All pass constant for mel-cepstrum analysis. "
+            "If None, automatically set to the best dimension for the sampling."
+        ),
+    )
+    parser.add_argument(
+        "--n_fft",
+        default=1024,
+        type=int,
+        help="The number of FFT points.",
+    )
+    parser.add_argument(
+        "--n_shift",
+        default=256,
+        type=int,
+        help="The number of shift points.",
+    )
+    parser.add_argument(
+        "--f0min",
+        default=40,
+        type=int,
+        help="Minimum f0 value.",
+    )
+    parser.add_argument(
+        "--f0max",
+        default=800,
+        type=int,
+        help="Maximum f0 value.",
+    )
+    parser.add_argument(
+        "--nj",
+        default=16,
+        type=int,
+        help="Number of parallel jobs.",
+    )
+    parser.add_argument(
+        "--verbose",
+        default=1,
+        type=int,
+        help="Verbosity level. Higher is more logging.",
+    )
+    return parser
+
+
+def main():
+    """Run log-F0 RMSE calculation in parallel."""
+    args = get_parser().parse_args()
+
+    # logging info
+    if args.verbose > 1:
+        logging.basicConfig(
+            level=logging.DEBUG,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    elif args.verbose > 0:
+        logging.basicConfig(
+            level=logging.INFO,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    else:
+        logging.basicConfig(
+            level=logging.WARN,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+        logging.warning("Skip DEBUG/INFO messages")
+
+    # find files
+    if os.path.isdir(args.gen_wavdir_or_wavscp):
+        gen_files = sorted(find_files(args.gen_wavdir_or_wavscp))
+    else:
+        with open(args.gen_wavdir_or_wavscp) as f:
+            gen_files = [line.strip().split(None, 1)[1] for line in f.readlines()]
+        if gen_files[0].endswith("|"):
+            raise ValueError("Not supported wav.scp format.")
+    if os.path.isdir(args.gt_wavdir_or_wavscp):
+        gt_files = sorted(find_files(args.gt_wavdir_or_wavscp))
+    else:
+        with open(args.gt_wavdir_or_wavscp) as f:
+            gt_files = [line.strip().split(None, 1)[1] for line in f.readlines()]
+        if gt_files[0].endswith("|"):
+            raise ValueError("Not supported wav.scp format.")
+
+    # Get and divide list
+    if len(gen_files) == 0:
+        raise FileNotFoundError("Not found any generated audio files.")
+    if len(gen_files) > len(gt_files):
+        raise ValueError(
+            "#groundtruth files are less than #generated files "
+            f"(#gen={len(gen_files)} vs. #gt={len(gt_files)}). "
+            "Please check the groundtruth directory."
+        )
+    logging.info("The number of utterances = %d" % len(gen_files))
+    file_lists = np.array_split(gen_files, args.nj)
+    file_lists = [f_list.tolist() for f_list in file_lists]
+
+    # multi processing
+    with mp.Manager() as manager:
+        log_f0_rmse_dict = manager.dict()
+        processes = []
+        # for f in file_lists:
+        #     calculate(f, gt_files, args, log_f0_rmse_dict)
+        for f in file_lists:
+            p = mp.Process(target=calculate, args=(f, gt_files, args, log_f0_rmse_dict))
+            p.start()
+            processes.append(p)
+
+        # wait for all process
+        for p in processes:
+            p.join()
+
+        # convert to standard list
+        log_f0_rmse_dict = dict(log_f0_rmse_dict)
+
+        # calculate statistics
+        mean_log_f0_rmse = np.mean(np.array([v for v in log_f0_rmse_dict.values()]))
+        std_log_f0_rmse = np.std(np.array([v for v in log_f0_rmse_dict.values()]))
+        logging.info(f"Average: {mean_log_f0_rmse:.4f} ± {std_log_f0_rmse:.4f}")
+
+    # write results
+    if args.outdir is None:
+        if os.path.isdir(args.gen_wavdir_or_wavscp):
+            args.outdir = args.gen_wavdir_or_wavscp
+        else:
+            args.outdir = os.path.dirname(args.gen_wavdir_or_wavscp)
+    os.makedirs(args.outdir, exist_ok=True)
+    with open(f"{args.outdir}/utt2log_f0_rmse", "w") as f:
+        for utt_id in sorted(log_f0_rmse_dict.keys()):
+            log_f0_rmse = log_f0_rmse_dict[utt_id]
+            f.write(f"{utt_id} {log_f0_rmse:.4f}\n")
+    with open(f"{args.outdir}/log_f0_rmse_avg_result.txt", "w") as f:
+        f.write(f"#utterances: {len(gen_files)}\n")
+        f.write(f"Average: {mean_log_f0_rmse:.4f} ± {std_log_f0_rmse:.4f}")
+
+    logging.info("Successfully finished log-F0 RMSE evaluation.")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/utils/evaluate_mcd.py b/utils/evaluate_mcd.py
new file mode 100644
index 00000000..8ff7a225
--- /dev/null
+++ b/utils/evaluate_mcd.py
@@ -0,0 +1,332 @@
+#!/usr/bin/env python3
+
+# Copyright 2020 Wen-Chin Huang and Tomoki Hayashi
+#  Apache 2.0  (http://www.apache.org/licenses/LICENSE-2.0)
+
+"""Evaluate MCD between generated and groundtruth audios with SPTK-based mcep."""
+
+import argparse
+import fnmatch
+import logging
+import multiprocessing as mp
+import os
+from typing import Dict, List, Tuple
+
+import librosa
+import resampy
+import numpy as np
+import pysptk
+import soundfile as sf
+from fastdtw import fastdtw
+from scipy import spatial
+
+
+def find_files(
+    root_dir: str, query: List[str] = ["*.flac", "*.wav"], include_root_dir: bool = True
+) -> List[str]:
+    """Find files recursively.
+
+    Args:
+        root_dir (str): Root root_dir to find.
+        query (List[str]): Query to find.
+        include_root_dir (bool): If False, root_dir name is not included.
+
+    Returns:
+        List[str]: List of found filenames.
+
+    """
+    files = []
+    for root, dirnames, filenames in os.walk(root_dir, followlinks=True):
+        for q in query:
+            for filename in fnmatch.filter(filenames, q):
+                files.append(os.path.join(root, filename))
+    if not include_root_dir:
+        files = [file_.replace(root_dir + "/", "") for file_ in files]
+
+    return files
+
+
+def sptk_extract(
+    x: np.ndarray,
+    fs: int,
+    n_fft: int = 512,
+    n_shift: int = 256,
+    mcep_dim: int = 25,
+    mcep_alpha: float = 0.41,
+    is_padding: bool = False,
+) -> np.ndarray:
+    """Extract SPTK-based mel-cepstrum.
+
+    Args:
+        x (ndarray): 1D waveform array.
+        fs (int): Sampling rate
+        n_fft (int): FFT length in point (default=512).
+        n_shift (int): Shift length in point (default=256).
+        mcep_dim (int): Dimension of mel-cepstrum (default=25).
+        mcep_alpha (float): All pass filter coefficient (default=0.41).
+        is_padding (bool): Whether to pad the end of signal (default=False).
+
+    Returns:
+        ndarray: Mel-cepstrum with the size (N, n_fft).
+
+    """
+    # perform padding
+    if is_padding:
+        n_pad = n_fft - (len(x) - n_fft) % n_shift
+        x = np.pad(x, (0, n_pad), "reflect")
+
+    # get number of frames
+    n_frame = (len(x) - n_fft) // n_shift + 1
+
+    # get window function
+    win = pysptk.sptk.hamming(n_fft)
+
+    # check mcep and alpha
+    if mcep_dim is None or mcep_alpha is None:
+        mcep_dim, mcep_alpha = _get_best_mcep_params(fs)
+
+    # calculate spectrogram
+    mcep = [
+        pysptk.mcep(
+            x[n_shift * i : n_shift * i + n_fft] * win,
+            mcep_dim,
+            mcep_alpha,
+            eps=1e-6,
+            etype=1,
+        )
+        for i in range(n_frame)
+    ]
+
+    return np.stack(mcep)
+
+
+def _get_basename(path: str) -> str:
+    return os.path.splitext(os.path.split(path)[-1])[0]
+
+
+def _get_best_mcep_params(fs: int) -> Tuple[int, float]:
+    if fs == 16000:
+        return 23, 0.42
+    elif fs == 22050:
+        return 34, 0.45
+    elif fs == 24000:
+        return 34, 0.46
+    elif fs == 44100:
+        return 39, 0.53
+    elif fs == 48000:
+        return 39, 0.55
+    else:
+        raise ValueError(f"Not found the setting for {fs}.")
+
+
+def calculate(
+    file_list: List[str],
+    gt_file_list: List[str],
+    args: argparse.Namespace,
+    mcd_dict: Dict,
+):
+    """Calculate MCD."""
+    for i, gen_path in enumerate(file_list):
+        corresponding_list = list(
+            filter(lambda gt_path: _get_basename(gt_path) in gen_path, gt_file_list)
+        )
+        assert len(corresponding_list) == 1
+        gt_path = corresponding_list[0]
+        gt_basename = _get_basename(gt_path)
+
+        # load wav file as int16
+        gen_x, gen_fs = sf.read(gen_path, dtype="int16")
+        gt_x, gt_fs = sf.read(gt_path, dtype="int16")
+
+        fs = gen_fs
+        if gen_fs != gt_fs:
+            gt_x = librosa.resample(gt_x.astype(np.float64), orig_sr=gt_fs, target_sr=gen_fs)
+
+        # extract ground truth and converted features
+        gen_mcep = sptk_extract(
+            x=gen_x,
+            fs=fs,
+            n_fft=args.n_fft,
+            n_shift=args.n_shift,
+            mcep_dim=args.mcep_dim,
+            mcep_alpha=args.mcep_alpha,
+        )
+        gt_mcep = sptk_extract(
+            x=gt_x,
+            fs=fs,
+            n_fft=args.n_fft,
+            n_shift=args.n_shift,
+            mcep_dim=args.mcep_dim,
+            mcep_alpha=args.mcep_alpha,
+        )
+
+        # DTW
+        _, path = fastdtw(gen_mcep, gt_mcep, dist=spatial.distance.euclidean)
+        twf = np.array(path).T
+        gen_mcep_dtw = gen_mcep[twf[0]]
+        gt_mcep_dtw = gt_mcep[twf[1]]
+
+        # MCD
+        diff2sum = np.sum((gen_mcep_dtw - gt_mcep_dtw) ** 2, 1)
+        mcd = np.mean(10.0 / np.log(10.0) * np.sqrt(2 * diff2sum), 0)
+        logging.info(f"{gt_basename} {mcd:.4f}")
+        mcd_dict[gt_basename] = mcd
+
+
+def get_parser() -> argparse.Namespace:
+    """Get argument parser."""
+    parser = argparse.ArgumentParser(description="Evaluate Mel-cepstrum distortion.")
+    parser.add_argument(
+        "gen_wavdir_or_wavscp",
+        type=str,
+        help="Path of directory or wav.scp for generated waveforms.",
+    )
+    parser.add_argument(
+        "gt_wavdir_or_wavscp",
+        type=str,
+        help="Path of directory or wav.scp for ground truth waveforms.",
+    )
+    parser.add_argument(
+        "--outdir",
+        type=str,
+        help="Path of directory to write the results.",
+    )
+
+    # analysis related
+    parser.add_argument(
+        "--mcep_dim",
+        default=None,
+        type=int,
+        help=(
+            "Dimension of mel cepstrum coefficients. "
+            "If None, automatically set to the best dimension for the sampling."
+        ),
+    )
+    parser.add_argument(
+        "--mcep_alpha",
+        default=None,
+        type=float,
+        help=(
+            "All pass constant for mel-cepstrum analysis. "
+            "If None, automatically set to the best dimension for the sampling."
+        ),
+    )
+    parser.add_argument(
+        "--n_fft",
+        default=1024,
+        type=int,
+        help="The number of FFT points.",
+    )
+    parser.add_argument(
+        "--n_shift",
+        default=256,
+        type=int,
+        help="The number of shift points.",
+    )
+    parser.add_argument(
+        "--nj",
+        default=16,
+        type=int,
+        help="Number of parallel jobs.",
+    )
+    parser.add_argument(
+        "--verbose",
+        default=1,
+        type=int,
+        help="Verbosity level. Higher is more logging.",
+    )
+    return parser
+
+
+def main():
+    """Run MCD calculation in parallel."""
+    args = get_parser().parse_args()
+
+    # logging info
+    if args.verbose > 1:
+        logging.basicConfig(
+            level=logging.DEBUG,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    elif args.verbose > 0:
+        logging.basicConfig(
+            level=logging.INFO,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    else:
+        logging.basicConfig(
+            level=logging.WARN,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+        logging.warning("Skip DEBUG/INFO messages")
+
+    # find files
+    if os.path.isdir(args.gen_wavdir_or_wavscp):
+        gen_files = sorted(find_files(args.gen_wavdir_or_wavscp))
+    else:
+        with open(args.gen_wavdir_or_wavscp) as f:
+            gen_files = [line.strip().split(None, 1)[1] for line in f.readlines()]
+        if gen_files[0].endswith("|"):
+            raise ValueError("Not supported wav.scp format.")
+    if os.path.isdir(args.gt_wavdir_or_wavscp):
+        gt_files = sorted(find_files(args.gt_wavdir_or_wavscp))
+    else:
+        with open(args.gt_wavdir_or_wavscp) as f:
+            gt_files = [line.strip().split(None, 1)[1] for line in f.readlines()]
+        if gt_files[0].endswith("|"):
+            raise ValueError("Not supported wav.scp format.")
+
+    # Get and divide list
+    if len(gen_files) == 0:
+        raise FileNotFoundError("Not found any generated audio files.")
+    if len(gen_files) > len(gt_files):
+        raise ValueError(
+            "#groundtruth files are less than #generated files "
+            f"(#gen={len(gen_files)} vs. #gt={len(gt_files)}). "
+            "Please check the groundtruth directory."
+        )
+    logging.info("The number of utterances = %d" % len(gen_files))
+    file_lists = np.array_split(gen_files, args.nj)
+    file_lists = [f_list.tolist() for f_list in file_lists]
+
+    # multi processing
+    with mp.Manager() as manager:
+        mcd_dict = manager.dict()
+        processes = []
+        for f in file_lists:
+            p = mp.Process(target=calculate, args=(f, gt_files, args, mcd_dict))
+            p.start()
+            processes.append(p)
+
+        # wait for all process
+        for p in processes:
+            p.join()
+
+        # convert to standard list
+        mcd_dict = dict(mcd_dict)
+
+        # calculate statistics
+        mean_mcd = np.mean(np.array([v for v in mcd_dict.values()]))
+        std_mcd = np.std(np.array([v for v in mcd_dict.values()]))
+        logging.info(f"Average: {mean_mcd:.4f} ± {std_mcd:.4f}")
+
+    # write results
+    if args.outdir is None:
+        if os.path.isdir(args.gen_wavdir_or_wavscp):
+            args.outdir = args.gen_wavdir_or_wavscp
+        else:
+            args.outdir = os.path.dirname(args.gen_wavdir_or_wavscp)
+    os.makedirs(args.outdir, exist_ok=True)
+    with open(f"{args.outdir}/utt2mcd", "w") as f:
+        for utt_id in sorted(mcd_dict.keys()):
+            mcd = mcd_dict[utt_id]
+            f.write(f"{utt_id} {mcd:.4f}\n")
+    with open(f"{args.outdir}/mcd_avg_result.txt", "w") as f:
+        f.write(f"#utterances: {len(gen_files)}\n")
+        f.write(f"Average: {mean_mcd:.4f} ± {std_mcd:.4f}")
+
+    logging.info("Successfully finished MCD evaluation.")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/utils/evaluate_semitone.py b/utils/evaluate_semitone.py
new file mode 100644
index 00000000..6c63c6c7
--- /dev/null
+++ b/utils/evaluate_semitone.py
@@ -0,0 +1,356 @@
+#!/usr/bin/env python3
+
+# Copyright 2021 Wen-Chin Huang and Tomoki Hayashi
+# Copyright 2022 Shuai Guo
+#  Apache 2.0  (http://www.apache.org/licenses/LICENSE-2.0)
+
+"""Evaluate semitone ACC between generated and groundtruth audios based on World."""
+
+import argparse
+import fnmatch
+import logging
+import multiprocessing as mp
+import os
+from math import log2, pow
+from typing import Dict, List, Tuple
+
+import librosa
+import numpy as np
+import pysptk
+import pyworld as pw
+import soundfile as sf
+from fastdtw import fastdtw
+from scipy import spatial
+
+
+def _Hz2Semitone(freq):
+    """_Hz2Semitone."""
+    A4 = 440
+    C0 = A4 * pow(2, -4.75)
+    name = ["C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B"]
+
+    if freq == 0:
+        return "Sil"  # silence
+    else:
+        h = round(12 * log2(freq / C0))
+        octave = h // 12
+        n = h % 12
+        return name[n] + "_" + str(octave)
+
+
+def find_files(
+    root_dir: str, query: List[str] = ["*.flac", "*.wav"], include_root_dir: bool = True
+) -> List[str]:
+    """Find files recursively.
+
+    Args:
+        root_dir (str): Root root_dir to find.
+        query (List[str]): Query to find.
+        include_root_dir (bool): If False, root_dir name is not included.
+
+    Returns:
+        List[str]: List of found filenames.
+
+    """
+    files = []
+    for root, dirnames, filenames in os.walk(root_dir, followlinks=True):
+        for q in query:
+            for filename in fnmatch.filter(filenames, q):
+                files.append(os.path.join(root, filename))
+    if not include_root_dir:
+        files = [file_.replace(root_dir + "/", "") for file_ in files]
+
+    return files
+
+
+def world_extract(
+    x: np.ndarray,
+    fs: int,
+    f0min: int = 40,
+    f0max: int = 800,
+    n_fft: int = 512,
+    n_shift: int = 256,
+    mcep_dim: int = 25,
+    mcep_alpha: float = 0.41,
+) -> np.ndarray:
+    """Extract World-based acoustic features.
+
+    Args:
+        x (ndarray): 1D waveform array.
+        fs (int): Minimum f0 value (default=40).
+        f0 (int): Maximum f0 value (default=800).
+        n_shift (int): Shift length in point (default=256).
+        n_fft (int): FFT length in point (default=512).
+        n_shift (int): Shift length in point (default=256).
+        mcep_dim (int): Dimension of mel-cepstrum (default=25).
+        mcep_alpha (float): All pass filter coefficient (default=0.41).
+
+    Returns:
+        ndarray: Mel-cepstrum with the size (N, n_fft).
+        ndarray: F0 sequence (N,).
+
+    """
+    # extract features
+    x = x.astype(np.float64)
+    f0, time_axis = pw.harvest(
+        x,
+        fs,
+        f0_floor=f0min,
+        f0_ceil=f0max,
+        frame_period=n_shift / fs * 1000,
+    )
+    sp = pw.cheaptrick(x, f0, time_axis, fs, fft_size=n_fft)
+    if mcep_dim is None or mcep_alpha is None:
+        mcep_dim, mcep_alpha = _get_best_mcep_params(fs)
+    mcep = pysptk.sp2mc(sp, mcep_dim, mcep_alpha)
+
+    return mcep, f0
+
+
+def _get_basename(path: str) -> str:
+    return os.path.splitext(os.path.split(path)[-1])[0]
+
+
+def _get_best_mcep_params(fs: int) -> Tuple[int, float]:
+    if fs == 16000:
+        return 23, 0.42
+    elif fs == 22050:
+        return 34, 0.45
+    elif fs == 24000:
+        return 34, 0.46
+    elif fs == 44100:
+        return 39, 0.53
+    elif fs == 48000:
+        return 39, 0.55
+    else:
+        raise ValueError(f"Not found the setting for {fs}.")
+
+
+def calculate(
+    file_list: List[str],
+    gt_file_list: List[str],
+    args: argparse.Namespace,
+    semitone_acc_dict: Dict[str, float],
+):
+    """Calculate semitone ACC."""
+    for i, gen_path in enumerate(file_list):
+        corresponding_list = list(
+            filter(lambda gt_path: _get_basename(gt_path) in gen_path, gt_file_list)
+        )
+        assert len(corresponding_list) == 1
+        gt_path = corresponding_list[0]
+        gt_basename = _get_basename(gt_path)
+
+        # load wav file as int16
+        gen_x, gen_fs = sf.read(gen_path, dtype="int16")
+        gt_x, gt_fs = sf.read(gt_path, dtype="int16")
+
+        fs = gen_fs
+        if gen_fs != gt_fs:
+            gt_x = librosa.resample(gt_x.astype(np.float), gt_fs, gen_fs)
+
+        # extract ground truth and converted features
+        gen_mcep, gen_f0 = world_extract(
+            x=gen_x,
+            fs=fs,
+            f0min=args.f0min,
+            f0max=args.f0max,
+            n_fft=args.n_fft,
+            n_shift=args.n_shift,
+            mcep_dim=args.mcep_dim,
+            mcep_alpha=args.mcep_alpha,
+        )
+        gt_mcep, gt_f0 = world_extract(
+            x=gt_x,
+            fs=fs,
+            f0min=args.f0min,
+            f0max=args.f0max,
+            n_fft=args.n_fft,
+            n_shift=args.n_shift,
+            mcep_dim=args.mcep_dim,
+            mcep_alpha=args.mcep_alpha,
+        )
+
+        # DTW
+        _, path = fastdtw(gen_mcep, gt_mcep, dist=spatial.distance.euclidean)
+        twf = np.array(path).T
+        gen_f0_dtw = gen_f0[twf[0]]
+        gt_f0_dtw = gt_f0[twf[1]]
+
+        # Semitone ACC
+        semitone_GT = np.array([_Hz2Semitone(_f0) for _f0 in gt_f0_dtw])
+        semitone_predict = np.array([_Hz2Semitone(_f0) for _f0 in gen_f0_dtw])
+        semitone_ACC = float((semitone_GT == semitone_predict).sum()) / len(semitone_GT)
+        semitone_acc_dict[gt_basename] = semitone_ACC
+
+
+def get_parser() -> argparse.Namespace:
+    """Get argument parser."""
+    parser = argparse.ArgumentParser(description="Evaluate Mel-cepstrum distortion.")
+    parser.add_argument(
+        "gen_wavdir_or_wavscp",
+        type=str,
+        help="Path of directory or wav.scp for generated waveforms.",
+    )
+    parser.add_argument(
+        "gt_wavdir_or_wavscp",
+        type=str,
+        help="Path of directory or wav.scp for ground truth waveforms.",
+    )
+    parser.add_argument(
+        "--outdir",
+        type=str,
+        help="Path of directory to write the results.",
+    )
+
+    # analysis related
+    parser.add_argument(
+        "--mcep_dim",
+        default=None,
+        type=int,
+        help=(
+            "Dimension of mel cepstrum coefficients. "
+            "If None, automatically set to the best dimension for the sampling."
+        ),
+    )
+    parser.add_argument(
+        "--mcep_alpha",
+        default=None,
+        type=float,
+        help=(
+            "All pass constant for mel-cepstrum analysis. "
+            "If None, automatically set to the best dimension for the sampling."
+        ),
+    )
+    parser.add_argument(
+        "--n_fft",
+        default=1024,
+        type=int,
+        help="The number of FFT points.",
+    )
+    parser.add_argument(
+        "--n_shift",
+        default=256,
+        type=int,
+        help="The number of shift points.",
+    )
+    parser.add_argument(
+        "--f0min",
+        default=40,
+        type=int,
+        help="Minimum f0 value.",
+    )
+    parser.add_argument(
+        "--f0max",
+        default=800,
+        type=int,
+        help="Maximum f0 value.",
+    )
+    parser.add_argument(
+        "--nj",
+        default=16,
+        type=int,
+        help="Number of parallel jobs.",
+    )
+    parser.add_argument(
+        "--verbose",
+        default=1,
+        type=int,
+        help="Verbosity level. Higher is more logging.",
+    )
+    return parser
+
+
+def main():
+    """Run semitone ACC calculation in parallel."""
+    args = get_parser().parse_args()
+
+    # logging info
+    if args.verbose > 1:
+        logging.basicConfig(
+            level=logging.DEBUG,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    elif args.verbose > 0:
+        logging.basicConfig(
+            level=logging.INFO,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    else:
+        logging.basicConfig(
+            level=logging.WARN,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+        logging.warning("Skip DEBUG/INFO messages")
+
+    # find files
+    if os.path.isdir(args.gen_wavdir_or_wavscp):
+        gen_files = sorted(find_files(args.gen_wavdir_or_wavscp))
+    else:
+        with open(args.gen_wavdir_or_wavscp) as f:
+            gen_files = [line.strip().split(None, 1)[1] for line in f.readlines()]
+        if gen_files[0].endswith("|"):
+            raise ValueError("Not supported wav.scp format.")
+    if os.path.isdir(args.gt_wavdir_or_wavscp):
+        gt_files = sorted(find_files(args.gt_wavdir_or_wavscp))
+    else:
+        with open(args.gt_wavdir_or_wavscp) as f:
+            gt_files = [line.strip().split(None, 1)[1] for line in f.readlines()]
+        if gt_files[0].endswith("|"):
+            raise ValueError("Not supported wav.scp format.")
+
+    # Get and divide list
+    if len(gen_files) == 0:
+        raise FileNotFoundError("Not found any generated audio files.")
+    if len(gen_files) > len(gt_files):
+        raise ValueError(
+            "#groundtruth files are less than #generated files "
+            f"(#gen={len(gen_files)} vs. #gt={len(gt_files)}). "
+            "Please check the groundtruth directory."
+        )
+    logging.info("The number of utterances = %d" % len(gen_files))
+    file_lists = np.array_split(gen_files, args.nj)
+    file_lists = [f_list.tolist() for f_list in file_lists]
+
+    # multi processing
+    with mp.Manager() as manager:
+        semitone_acc_dict = manager.dict()
+        processes = []
+        for f in file_lists:
+            p = mp.Process(
+                target=calculate, args=(f, gt_files, args, semitone_acc_dict)
+            )
+            p.start()
+            processes.append(p)
+
+        # wait for all process
+        for p in processes:
+            p.join()
+
+        # convert to standard list
+        semitone_acc_dict = dict(semitone_acc_dict)
+
+        # calculate statistics
+        mean_semitone_acc = np.mean(np.array([v for v in semitone_acc_dict.values()]))
+        logging.info(f"Average - Semitone_ACC: {mean_semitone_acc*100:.2f}%")
+
+    # write results
+    if args.outdir is None:
+        if os.path.isdir(args.gen_wavdir_or_wavscp):
+            args.outdir = args.gen_wavdir_or_wavscp
+        else:
+            args.outdir = os.path.dirname(args.gen_wavdir_or_wavscp)
+    os.makedirs(args.outdir, exist_ok=True)
+    with open(f"{args.outdir}/utt2semitone_acc", "w") as f:
+        for utt_id in sorted(semitone_acc_dict.keys()):
+            semitone_ACC = semitone_acc_dict[utt_id]
+            f.write(f"{utt_id} {semitone_ACC*100:.2f}%\n")
+    with open(f"{args.outdir}/semitone_acc_avg_result.txt", "w") as f:
+        f.write(f"#utterances: {len(gen_files)}\n")
+        f.write(f"Average: {mean_semitone_acc*100:.2f}%")
+
+    logging.info("Successfully finished semitone ACC evaluation.")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/utils/evaluate_vuv.py b/utils/evaluate_vuv.py
new file mode 100644
index 00000000..0e3499c9
--- /dev/null
+++ b/utils/evaluate_vuv.py
@@ -0,0 +1,345 @@
+#!/usr/bin/env python3
+
+# Copyright 2021 Wen-Chin Huang and Tomoki Hayashi
+# Copyright 2022 Shuai Guo
+#  Apache 2.0  (http://www.apache.org/licenses/LICENSE-2.0)
+
+"""Evaluate VUV error between generated and groundtruth audios based on World."""
+
+import argparse
+import fnmatch
+import logging
+import multiprocessing as mp
+import os
+from typing import Dict, List, Tuple
+
+import librosa
+import numpy as np
+import pysptk
+import pyworld as pw
+import soundfile as sf
+from fastdtw import fastdtw
+from scipy import spatial
+
+
+def _Hz2Flag(freq):
+    if freq == 0:
+        return False
+    else:
+        return True
+
+
+def find_files(
+    root_dir: str, query: List[str] = ["*.flac", "*.wav"], include_root_dir: bool = True
+) -> List[str]:
+    """Find files recursively.
+
+    Args:
+        root_dir (str): Root root_dir to find.
+        query (List[str]): Query to find.
+        include_root_dir (bool): If False, root_dir name is not included.
+
+    Returns:
+        List[str]: List of found filenames.
+
+    """
+    files = []
+    for root, dirnames, filenames in os.walk(root_dir, followlinks=True):
+        for q in query:
+            for filename in fnmatch.filter(filenames, q):
+                files.append(os.path.join(root, filename))
+    if not include_root_dir:
+        files = [file_.replace(root_dir + "/", "") for file_ in files]
+
+    return files
+
+
+def world_extract(
+    x: np.ndarray,
+    fs: int,
+    f0min: int = 40,
+    f0max: int = 800,
+    n_fft: int = 512,
+    n_shift: int = 256,
+    mcep_dim: int = 25,
+    mcep_alpha: float = 0.41,
+) -> np.ndarray:
+    """Extract World-based acoustic features.
+
+    Args:
+        x (ndarray): 1D waveform array.
+        fs (int): Minimum f0 value (default=40).
+        f0 (int): Maximum f0 value (default=800).
+        n_shift (int): Shift length in point (default=256).
+        n_fft (int): FFT length in point (default=512).
+        n_shift (int): Shift length in point (default=256).
+        mcep_dim (int): Dimension of mel-cepstrum (default=25).
+        mcep_alpha (float): All pass filter coefficient (default=0.41).
+
+    Returns:
+        ndarray: Mel-cepstrum with the size (N, n_fft).
+        ndarray: F0 sequence (N,).
+
+    """
+    # extract features
+    x = x.astype(np.float64)
+    f0, time_axis = pw.harvest(
+        x,
+        fs,
+        f0_floor=f0min,
+        f0_ceil=f0max,
+        frame_period=n_shift / fs * 1000,
+    )
+    sp = pw.cheaptrick(x, f0, time_axis, fs, fft_size=n_fft)
+    if mcep_dim is None or mcep_alpha is None:
+        mcep_dim, mcep_alpha = _get_best_mcep_params(fs)
+    mcep = pysptk.sp2mc(sp, mcep_dim, mcep_alpha)
+
+    return mcep, f0
+
+
+def _get_basename(path: str) -> str:
+    return os.path.splitext(os.path.split(path)[-1])[0]
+
+
+def _get_best_mcep_params(fs: int) -> Tuple[int, float]:
+    if fs == 16000:
+        return 23, 0.42
+    elif fs == 22050:
+        return 34, 0.45
+    elif fs == 24000:
+        return 34, 0.46
+    elif fs == 44100:
+        return 39, 0.53
+    elif fs == 48000:
+        return 39, 0.55
+    else:
+        raise ValueError(f"Not found the setting for {fs}.")
+
+
+def calculate(
+    file_list: List[str],
+    gt_file_list: List[str],
+    args: argparse.Namespace,
+    vuv_err_dict: Dict[str, float],
+):
+    """Calculate VUV error."""
+    for i, gen_path in enumerate(file_list):
+        corresponding_list = list(
+            filter(lambda gt_path: _get_basename(gt_path) in gen_path, gt_file_list)
+        )
+        assert len(corresponding_list) == 1
+        gt_path = corresponding_list[0]
+        gt_basename = _get_basename(gt_path)
+
+        # load wav file as int16
+        gen_x, gen_fs = sf.read(gen_path, dtype="int16")
+        gt_x, gt_fs = sf.read(gt_path, dtype="int16")
+
+        fs = gen_fs
+        if gen_fs != gt_fs:
+            gt_x = librosa.resample(gt_x.astype(np.float64), orig_sr=gt_fs, target_sr=gen_fs)
+
+        # extract ground truth and converted features
+        gen_mcep, gen_f0 = world_extract(
+            x=gen_x,
+            fs=fs,
+            f0min=args.f0min,
+            f0max=args.f0max,
+            n_fft=args.n_fft,
+            n_shift=args.n_shift,
+            mcep_dim=args.mcep_dim,
+            mcep_alpha=args.mcep_alpha,
+        )
+        gt_mcep, gt_f0 = world_extract(
+            x=gt_x,
+            fs=fs,
+            f0min=args.f0min,
+            f0max=args.f0max,
+            n_fft=args.n_fft,
+            n_shift=args.n_shift,
+            mcep_dim=args.mcep_dim,
+            mcep_alpha=args.mcep_alpha,
+        )
+
+        # DTW
+        _, path = fastdtw(gen_mcep, gt_mcep, dist=spatial.distance.euclidean)
+        twf = np.array(path).T
+        gen_f0_dtw = gen_f0[twf[0]]
+        gt_f0_dtw = gt_f0[twf[1]]
+
+        # VUV ERR
+        vuv_GT = np.array([_Hz2Flag(_f0) for _f0 in gt_f0_dtw])
+        vuv_predict = np.array([_Hz2Flag(_f0) for _f0 in gen_f0_dtw])
+        vuv_ERR = float((vuv_GT != vuv_predict).sum()) / len(vuv_GT)
+        vuv_err_dict[gt_basename] = vuv_ERR
+
+
+def get_parser() -> argparse.Namespace:
+    """Get argument parser."""
+    parser = argparse.ArgumentParser(description="Evaluate Mel-cepstrum distortion.")
+    parser.add_argument(
+        "gen_wavdir_or_wavscp",
+        type=str,
+        help="Path of directory or wav.scp for generated waveforms.",
+    )
+    parser.add_argument(
+        "gt_wavdir_or_wavscp",
+        type=str,
+        help="Path of directory or wav.scp for ground truth waveforms.",
+    )
+    parser.add_argument(
+        "--outdir",
+        type=str,
+        help="Path of directory to write the results.",
+    )
+
+    # analysis related
+    parser.add_argument(
+        "--mcep_dim",
+        default=None,
+        type=int,
+        help=(
+            "Dimension of mel cepstrum coefficients. "
+            "If None, automatically set to the best dimension for the sampling."
+        ),
+    )
+    parser.add_argument(
+        "--mcep_alpha",
+        default=None,
+        type=float,
+        help=(
+            "All pass constant for mel-cepstrum analysis. "
+            "If None, automatically set to the best dimension for the sampling."
+        ),
+    )
+    parser.add_argument(
+        "--n_fft",
+        default=1024,
+        type=int,
+        help="The number of FFT points.",
+    )
+    parser.add_argument(
+        "--n_shift",
+        default=256,
+        type=int,
+        help="The number of shift points.",
+    )
+    parser.add_argument(
+        "--f0min",
+        default=40,
+        type=int,
+        help="Minimum f0 value.",
+    )
+    parser.add_argument(
+        "--f0max",
+        default=800,
+        type=int,
+        help="Maximum f0 value.",
+    )
+    parser.add_argument(
+        "--nj",
+        default=16,
+        type=int,
+        help="Number of parallel jobs.",
+    )
+    parser.add_argument(
+        "--verbose",
+        default=1,
+        type=int,
+        help="Verbosity level. Higher is more logging.",
+    )
+    return parser
+
+
+def main():
+    """Run VUV error calculation in parallel."""
+    args = get_parser().parse_args()
+
+    # logging info
+    if args.verbose > 1:
+        logging.basicConfig(
+            level=logging.DEBUG,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    elif args.verbose > 0:
+        logging.basicConfig(
+            level=logging.INFO,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+    else:
+        logging.basicConfig(
+            level=logging.WARN,
+            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
+        )
+        logging.warning("Skip DEBUG/INFO messages")
+
+    # find files
+    if os.path.isdir(args.gen_wavdir_or_wavscp):
+        gen_files = sorted(find_files(args.gen_wavdir_or_wavscp))
+    else:
+        with open(args.gen_wavdir_or_wavscp) as f:
+            gen_files = [line.strip().split(None, 1)[1] for line in f.readlines()]
+        if gen_files[0].endswith("|"):
+            raise ValueError("Not supported wav.scp format.")
+    if os.path.isdir(args.gt_wavdir_or_wavscp):
+        gt_files = sorted(find_files(args.gt_wavdir_or_wavscp))
+    else:
+        with open(args.gt_wavdir_or_wavscp) as f:
+            gt_files = [line.strip().split(None, 1)[1] for line in f.readlines()]
+        if gt_files[0].endswith("|"):
+            raise ValueError("Not supported wav.scp format.")
+
+    # Get and divide list
+    if len(gen_files) == 0:
+        raise FileNotFoundError("Not found any generated audio files.")
+    if len(gen_files) > len(gt_files):
+        raise ValueError(
+            "#groundtruth files are less than #generated files "
+            f"(#gen={len(gen_files)} vs. #gt={len(gt_files)}). "
+            "Please check the groundtruth directory."
+        )
+    logging.info("The number of utterances = %d" % len(gen_files))
+    file_lists = np.array_split(gen_files, args.nj)
+    file_lists = [f_list.tolist() for f_list in file_lists]
+
+    # multi processing
+    with mp.Manager() as manager:
+        vuv_err_dict = manager.dict()
+        processes = []
+        for f in file_lists:
+            p = mp.Process(target=calculate, args=(f, gt_files, args, vuv_err_dict))
+            p.start()
+            processes.append(p)
+
+        # wait for all process
+        for p in processes:
+            p.join()
+
+        # convert to standard list
+        vuv_err_dict = dict(vuv_err_dict)
+
+        # calculate statistics
+        mean_vuv_err = np.mean(np.array([v for v in vuv_err_dict.values()]))
+        logging.info(f"Average - VUV_ERROR: {mean_vuv_err*100:.2f}%")
+
+    # write results
+    if args.outdir is None:
+        if os.path.isdir(args.gen_wavdir_or_wavscp):
+            args.outdir = args.gen_wavdir_or_wavscp
+        else:
+            args.outdir = os.path.dirname(args.gen_wavdir_or_wavscp)
+    os.makedirs(args.outdir, exist_ok=True)
+    with open(f"{args.outdir}/utt2vuv_error", "w") as f:
+        for utt_id in sorted(vuv_err_dict.keys()):
+            vuv_ERR = vuv_err_dict[utt_id]
+            f.write(f"{utt_id} {vuv_ERR*100:.2f}%\n")
+    with open(f"{args.outdir}/vuv_error_avg_result.txt", "w") as f:
+        f.write(f"#utterances: {len(gen_files)}\n")
+        f.write(f"Average: {mean_vuv_err*100:.2f}%")
+
+    logging.info("Successfully finished VUV error evaluation.")
+
+
+if __name__ == "__main__":
+    main()