# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Using Spectrogram to visualize the nonstationary frequency response of the signal
#
# ## 0. Introduction
# Given signal of length L, and the STFT parameters:
# 1. Window length, $M$
# 2. Shift/stride, $R$ ($ 1 \le R \le M $, for no loss of information)
# 3. FFT size, $N$ ($N\ge M$, for our purpose, $ N=M $)
# Then # segments, $K$, will be $ K=\lfloor (L-M)/R \rfloor$
# In our problem, our data samples have $L=128$, which limits our options for window length. If we choose a large $M$ necessary for finer resolution of the frequency components (say $M\ge L$ with zero-padding or over-sampling), we would lose the temporal information of when the frequency peaks occur. So we will make the following tradeoff: $N=M=32$, $R=2$. Furthermore prefix and suffix $M/2$ samples to the signal frame to also fully incorporate the spectral behavior at the edges (when using tapered windows), thus $L'=128+N=160$. With these parameters and adjustments, $K=64$. Thus the inputs to our CNN classifier will be $(M/2+1)\times K=17\times 64$ spectrogram images.

import os
import numpy as np
from scipy import signal
from scipy.fftpack import fft
from scipy.fft import fftshift
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.patches as patches
import itertools


# ##. Loading the UCI HAR dataset into an numpy ndarray
# Download dataset from https://archive.ics.uci.edu/ml/datasets/Human+Activity+Recognition+Using+Smartphones

# +
def read_signals(filename):
    with open(filename, 'r') as fp:
        data = fp.read().splitlines()
        data = map(lambda x: x.strip().split(), data)
        data = [list(map(float, line)) for line in data]
    return data

DATA_FOLDER = '../datasets/UCI HAR Dataset/'
INPUT_FOLDER_TRAIN = DATA_FOLDER+'train/Inertial Signals/'
INPUT_FILES_TRAIN = ['body_acc_x_train.txt']

train_signals = []

for input_file in INPUT_FILES_TRAIN:
    sig = read_signals(INPUT_FOLDER_TRAIN + input_file)
    train_signals.append(sig)
train_signals = np.transpose(train_signals, (1, 2, 0))
train_signals = np.squeeze(train_signals)
[no_signals_train, N] = np.shape(train_signals)
time = np.arange(N)
xsignal = train_signals[10]


# -

# ## 1. Plot one signal from the UCI HAR dataset

# +
# First lets plot a signal frame together with its time-average

def get_ave_values(xvalues, yvalues, n = 5):
    signal_length = len(xvalues)
    if signal_length % n == 0:
        padding_length = 0
    else:
        padding_length = n - signal_length//n % n
    xarr = np.array(xvalues)
    yarr = np.array(yvalues)
    xarr.resize(signal_length//n, n)
    yarr.resize(signal_length//n, n)
    xarr_reshaped = xarr.reshape((-1,n))
    yarr_reshaped = yarr.reshape((-1,n))
    x_ave = xarr_reshaped[:,0]
    y_ave = np.nanmean(yarr_reshaped, axis=1)
    return x_ave, y_ave

def plot_signal_plus_average(ax, time, sig, average_over = 5):
    time_ave, signal_ave = get_ave_values(time, sig, average_over)
    ax.plot(time, sig, label='signal')
    ax.plot(time_ave, signal_ave, label = 'time average (n={})'.format(5))
    ax.set_xlim([time[0], time[-1]])
    ax.set_ylabel('Amplitude', fontsize=16)
    ax.set_title('Signal + Time Average', fontsize=16)
    ax.legend(loc='upper right')

fig, ax = plt.subplots(figsize=(12,3))
plot_signal_plus_average(ax, time, xsignal, average_over = 3)
plt.show()


# -

# ## 2. Plot the Fourier Transform of the body acceleration signal

# +
def get_fft_values(y_values, T, N, f_s):
    N2 = 2 ** (int(np.log2(N)) + 1) # round up to next highest power of 2
    f_values = np.linspace(0.0, 1.0/(2.0*T), N2//2)
    fft_values_ = fft(y_values)
    fft_values = 2.0/N2 * np.abs(fft_values_[0:N2//2])
    return f_values, fft_values

def plot_fft_plus_power(ax, time, sig, yticks=None, ylim=None):
    dt = time[1] - time[0]
    N = len(sig)
    fs = 1/dt
    
    variance = np.std(sig)**2
    f_values, fft_values = get_fft_values(sig, dt, N, fs)
    fft_power = variance * abs(fft_values) ** 2
    ax.plot(f_values, fft_values, 'r-', label='Fourier Transform')
    ax.plot(f_values, fft_power, 'k--', linewidth=1, label='FFT Power Spectrum')
    ax.legend()

fig, ax = plt.subplots(figsize=(12,3))
ax.set_xlabel('Frequency [Hz / year]', fontsize=18)
ax.set_ylabel('Amplitude', fontsize=18)
plot_fft_plus_power(ax, time, xsignal)
plt.show()


# -

# ## 3. Plot the Spectrogram

# +
def plot_spectrogram(ax, sig, M, noverlap, fs=1.0, windowname = ('hamming'),
                 cmap = plt.cm.seismic, title = '', ylabel = '', xlabel = ''):
    # get the window taps
    win = signal.get_window(windowname,M,False)
    # prefix/suffix
    pref = sig[-int(M/2):]*win[0:int(M/2)]
    suf = sig[0:int(M/2)]*win[-int(M/2):]

    sig = np.concatenate((pref, sig))
    sig = np.concatenate((sig,suf))
    f, t, Sxx = signal.spectrogram(sig, fs, window=win, nperseg=M, noverlap=noverlap)
    ax.pcolormesh(t, f, Sxx, shading='nearest',cmap=cmap)
    ax.set_title(title, fontsize=14)
    ax.set_ylabel(ylabel, fontsize=12)
    ax.set_xlabel(xlabel, fontsize=12)
    #Sxx = Sxx[:,0:-1]
    return f,t,Sxx

fig = plt.figure(figsize=(12,12))
spec = gridspec.GridSpec(ncols=6, nrows=8)
top_left_ax = fig.add_subplot(spec[0:3, 0:3])
top_right_ax = fig.add_subplot(spec[0:3, 3:6])
bottom_ax = fig.add_subplot(spec[4:, 1:5])

ylabel = 'Frequency'
xlabel = 'Time'
M = 32
noverlap = M-2
title = 'Power Spectrum (Hann)'
cmap = cmap=plt.get_cmap('jet')
f,t,Sxx1 = plot_spectrogram(top_left_ax, xsignal, M, noverlap, 
                           windowname='hann', cmap=cmap, xlabel=xlabel, ylabel=ylabel, title=title)
title = 'Power Spectrum (Boxcar)'
f,t,Sxx2 = plot_spectrogram(top_right_ax, xsignal, M, noverlap, 
                           windowname='boxcar', cmap=cmap, xlabel=xlabel, ylabel=ylabel, title=title)

#f,t,Sxx = plot_spectrogram(top_left_ax, xsignal, M, noverlap, xlabel=xlabel, ylabel=ylabel, title=title)
title = 'min(boxcar,hann)'
bottom_ax.pcolormesh(t, f, np.minimum(Sxx1,Sxx2), shading='nearest',cmap=cmap)
bottom_ax.set_title(title, fontsize=14)
bottom_ax.set_ylabel(ylabel, fontsize=12)
bottom_ax.set_xlabel(xlabel, fontsize=12)

plt.show()
# -



StarcoderdataJupyterScriptsDedupFiltered