voxelwiseencoding/make_stim_spectrogram_plots.py at master · ANCPLabOldenburg/voxelwiseencoding · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
# -*- coding: utf-8 -*-

#!/usr/bin/env python
import numpy as np
import glob
import os.path
#import json
import re
#import scipy.io.wavfile as wav
from librosa import load
from scipy.stats import gamma
from scipy.signal import spectrogram, welch, hilbert, fftconvolve, decimate
import matplotlib.pyplot as plt

def make_hrf(time_points, overshoot_delay=6, undershoot_delay=12,
             amplitiude_factor_undershoot=0.35):
    """ Make a simple hrf for the time points passed. """
    hrf = gamma.pdf(time_points, overshoot_delay) - \
        gamma.pdf(time_points, undershoot_delay) * amplitiude_factor_undershoot
    # Normalize to integral 1
    return hrf / sum(hrf)

def plot_timeseries(t, x, title, save):
    plt.plot(t, x)
    plt.xlabel('Time (s)')
    plt.title(title)
    plt.gcf().tight_layout()
#    plt.show()
    plt.savefig(save)
    plt.close()

if __name__ == "__main__":
    output_folder = '/data2/azubaidi/ForrestGumpHearingLoss/BIDS_ForrGump/derivatives/fmriprep/'
    output_sorted = True
    use_noise = False

    if use_noise:
        stim_folder = '/data2/azubaidi/ForrestGumpHearingLoss/BIDS_ForrGump/'\
                      +'sourcedata/stimuli/RecordedScannerNoise/'
    else:
        stim_folder = '/data2/azubaidi/ForrestGumpHearingLoss/BIDS_ForrGump/'\
                      +'sourcedata/stimuli/PresentedStimuli/'
#    stim_extension_old = 'recstimuli'
#    stim_extension = 'stim'
#    recording_extension = 'recording-rec_'
#    recording_extension = ''
#    subjects = ["01","02","03","04","05","06","07","08","09","10"]
    subjects = ["03"]
#    subjects = ["09"]
#    sessions = ["01","02","03"]
    sessions = ["01"]
    # Name of the folders where unincluded runs (01 and 08) go for each session
    unincluded_runs = {"01":"1st","02":"2nd","03":"3rd"}
    for subject in subjects:
        for session in sessions:
            subses_stim_folder = os.path.join(stim_folder, f"/ses-{session}/")
#            subses_stim_folder = os.path.join(stim_folder, f"sub-{subject}/ses-{session}/")
            wav_files = glob.glob(subses_stim_folder + "*.wav")
            if output_sorted:
                subses_output_folder_sorted = os.path.join(output_folder, f"sub-{subject}/ses-%s/func/")
            else:
                subses_output_folder = os.path.join(output_folder, f"sub-{subject}/ses-{session}/func/")

            for wav_file in wav_files:
#                if 'run-01' in wav_file or 'run-08' in wav_file:
##                    print("Skipping ",wav_file)
#                    continue

                outfile_base = os.path.basename(wav_file).split('.')[0]
                if use_noise:
                    outfile_base += '_stim'
#                outfile_base = outfile_base.replace(stim_extension_old,
#                                                    recording_extension+stim_extension)
                if output_sorted:
                    if 'run-01' in wav_file or 'run-08' in wav_file:
                        condition = unincluded_runs[session]
                    elif 'acq-CS' in wav_file:
                        condition = 'CS'
                    elif 'acq-N4' in wav_file:
                        condition = 'N4'
                    elif 'acq-S2' in wav_file:
                        condition = 'S2'
                    else:
                        raise Exception('No valid condition tag found in wav '
                                        +'file name. Valid tags: CS, N4, S2')
                    subses_output_folder = subses_output_folder_sorted % condition
                    outfile_base = re.sub(r'ses-[0-9]*','ses-'+condition,outfile_base)
                output_dir = os.path.join(subses_output_folder, "stim/")
                if not os.path.exists(output_dir):
                    os.mkdir(output_dir)
                output_dir += outfile_base
                # Found a memory leak in scipy.io.wavfile when reading a 16bit stereo
                # wav file: sys.getsizeof and numpys .__sizeof__ return 112
                # when the actual size is 159293432
                # garbage collection seems to fail and parts of the array are
                # overwritten while in use
                print("Converting", wav_file)
                sig2, Fs = load(wav_file, sr=None, mono=False, dtype=np.float64)
#                Fs, sig = wav.read(wav_file)
                # For some reason sig only has size 96 when it should have
                # size 155771172. Recasting the array as a different type
                # or performing an operation like mean seems to fix it.
                if sig2.ndim > 1:
                    sig = np.mean(sig2,axis=0,dtype=np.float32) # convert a WAV from stereo to mono
                else:
                    sig = sig2.copy()
                t = np.arange(0, len(sig)) / Fs

#                plot_timeseries(t, sig, "Wavfile timeseries", output_dir+"_timeseries.png")
                envelope = np.abs(hilbert(sig))
#                plot_timeseries(t, envelope, "Envelope", output_dir+"_envelope.png")

                TR = 0.85
                f, t, Sxx = spectrogram(envelope, fs=1/TR, nperseg=64)
                plt.pcolormesh(t, f, Sxx, shading='gourad')
                plt.xlabel('Time (s)')
                plt.ylabel('Frequency (Hz)')
                plt.title('Envelope spectrogram')
                plt.savefig(output_dir+'_envelope_spectrogram_nohrf.png')
#                plt.show()
                plt.close()

                # Make hemodynamic response function
#                sample_duration = 1 / Fs
#                time_points = np.arange(0, 25, sample_duration) # 25 seconds long
#                hrf = make_hrf(time_points)
#                plot_timeseries(time_points, hrf, "HRF", output_dir+"_hrf.png")
                # Convolve envelope with HRF
#                envelope_convolved = fftconvolve(envelope, hrf, mode='valid')
#                t = np.arange(0, len(envelope_convolved)) / Fs
#                plot_timeseries(t, envelope_convolved, "Envelope convolved with HRF",
#                                output_dir+"_envelope_convolved.png")
#                # Downsampling from wav file sampling rate 44100 Hz to
#                # fMRI sampling rate 1/TR (~1.176 Hz).
#                # Decimation factor: 44100 * TR = 37485 = 3^2 * 5 * 7^2 * 17
#                envelope_convolved = decimate(envelope_convolved,17)
#                envelope_convolved = decimate(envelope_convolved,9)
#                envelope_convolved = decimate(envelope_convolved,7)
#                envelope_convolved = decimate(envelope_convolved,7)
#                envelope_convolved = decimate(envelope_convolved,5)
#                # For some reason envelope_convolved only has size 96 when it
#                # should have size 8176.
#                envelope_convolved = np.asarray(envelope_convolved,dtype=np.float32)
#                t = np.arange(0, len(envelope_convolved)) * TR
#                plot_timeseries(t, envelope_convolved,
#                                "Envelope convolved with HRF and downsampled to fMRI Fs",
#                                output_dir+"_envelope_convolved_decimated.png")
#
#                f, t, Sxx = spectrogram(envelope_convolved, fs=1/TR, nperseg=64)
#                plt.pcolormesh(t, f, Sxx, shading='gourad')
#                plt.xlabel('Time (s)')
#                plt.ylabel('Frequency (Hz)')
#                plt.title('Envelope spectrogram')
#                plt.savefig(output_dir+'_envelope_spectrogram.png')
##                plt.show()
#                plt.close()
#
#                f, Pxx = welch(envelope_convolved, fs=1/TR)
#                plt.semilogy(f, Pxx)
#                plt.xlabel('Frequency (Hz)')
##                plt.ylabel('dB')
#                plt.title('Envelope Welch PSD')
#                plt.savefig(output_dir+'_envelope_periodogram_welch.png')
##                plt.show()
#                plt.close()
##                break
#                del sig, envelope, envelope_convolved
#                import gc
#                print(gc.collect())