FDD-task

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from mpl_toolkits.mplot3d import Axes3D
from tkinter import Tk, filedialog

# --- CONFIGURATION ---
SAMPLING_FREQUENCY = 250.0 
N_PER_SEG = 8192  # High resolution to match Artemis

# --- HELPER FUNCTIONS ---

def load_data():
    print("Select your Excel/CSV file...")
    root = Tk()
    root.withdraw() 
    root.attributes('-topmost', True)
    filename = filedialog.askopenfilename(title="Select Data File", filetypes=[("Data", "*.xlsx *.xls *.csv *.txt")])
    root.destroy()
    
    if not filename: return None
    
    try:
        if filename.endswith(('.xlsx', '.xls')):
            df = pd.read_excel(filename, header=0)
        else:
            try:
                df = pd.read_csv(filename, sep=None, engine='python', header=0)
            except:
                df = pd.read_csv(filename, sep=r'\s+', header=None)
                
        # Ensure numeric
        df = df.apply(pd.to_numeric, errors='coerce').dropna()
        print(f"Data Loaded: {df.shape}")
        return df.values
    except Exception as e:
        print(f"Error: {e}")
        return None

def compute_fdd(data, fs, nperseg):
    # Preprocess: Detrend and Bandpass (0.5 - 30 Hz)
    data = signal.detrend(data, axis=0, type='constant')
    sos = signal.butter(4, [0.5, 30], 'bandpass', fs=fs, output='sos')
    data = signal.sosfiltfilt(sos, data, axis=0)
    
    # Frequency Vector
    f, _ = signal.welch(data[:, 0], fs=fs, nperseg=nperseg)
    n_freqs = len(f)
    n_channels = data.shape[1]
    
    # CPSD Matrix
    G = np.zeros((n_freqs, n_channels, n_channels), dtype=complex)
    
    print("Calculating Spectral Matrices...")
    for i in range(n_channels):
        for j in range(i, n_channels):
            _, Pxy = signal.csd(data[:, i], data[:, j], fs=fs, nperseg=nperseg)
            G[:, i, j] = Pxy
            if i != j: G[:, j, i] = np.conj(Pxy)

    print("Performing SVD...")
    # Store ALL singular values
    S_all = np.zeros((n_freqs, n_channels))
    U = np.zeros((n_freqs, n_channels, n_channels), dtype=complex)
    
    for k in range(n_freqs):
        u, s, vh = np.linalg.svd(G[k, :, :])
        S_all[k, :] = s  # s is a vector of singular values (s1, s2, s3...)
        U[k, :, :] = u
        
    return f, S_all, U

def get_smart_peaks(freqs, s1):
    """
    Finds peaks using specific search zones to capture both strong and weak modes.
    """
    peaks, _ = signal.find_peaks(s1)
    
    candidates = []
    for p in peaks:
        f = freqs[p]
        if 1.5 < f < 20.0: # General bridge range
            candidates.append(p)
            
    # Define "Zones" where we expect modes (based on your Artemis data)
    zones = [
        (2.0, 3.0),   # Zone 1
        (3.0, 4.2),   # Zone 2
        (4.2, 5.5),   # Zone 3
        (6.0, 7.5),   # Zone 4
        (8.0, 10.0),  # Zone 5 (The Weak Mode)
        (10.5, 13.0)  # Zone 6 (The Strong Mode)
    ]
    
    final_peaks = []
    
    for (low, high) in zones:
        # Find candidates in this zone
        zone_candidates = [p for p in candidates if low <= freqs[p] <= high]
        
        if zone_candidates:
            # Pick the highest point in this specific zone
            best_p = max(zone_candidates, key=lambda x: s1[x])
            final_peaks.append(best_p)
            
    return np.array(final_peaks)

def plot_stick_model(u_matrix, peak_indices, freqs):
    """Plots 3D Stick Models."""
    fig = plt.figure(figsize=(15, 8))
    
    # Layout for 6 sensors (Linear Bridge Model)
    x_undef = np.linspace(0, 50, 6)
    y_undef = np.zeros(6)
    z_undef = np.zeros(6)
    
    for i, peak_idx in enumerate(peak_indices):
        if i >= 6: break 
        
        ax = fig.add_subplot(2, 3, i+1, projection='3d')
        
        freq = freqs[peak_idx]
        mode_shape = np.real(u_matrix[peak_idx, :, 0])
        
        # Normalize scale
        scale = 8.0 / np.max(np.abs(mode_shape))
        
        dx = mode_shape[0::3] * scale
        dy = mode_shape[1::3] * scale
        dz = mode_shape[2::3] * scale
        
        x_def = x_undef + dx
        y_def = y_undef + dy
        z_def = z_undef + dz
        
        # Plot Undeformed
        ax.plot(x_undef, y_undef, z_undef, 'k--', alpha=0.3)
        ax.scatter(x_undef, y_undef, z_undef, c='k', s=10, alpha=0.3)
        
        # Plot Deformed
        ax.plot(x_def, y_def, z_def, 'b-', linewidth=2)
        ax.scatter(x_def, y_def, z_def, c='r', s=30)
        
        # Connecting lines
        for j in range(6):
            ax.plot([x_undef[j], x_def[j]], [y_undef[j], y_def[j]], [z_undef[j], z_def[j]], 'r-', alpha=0.3)
            
        ax.set_title(f"Mode {i+1}: {freq:.3f} Hz", fontsize=10, fontweight='bold')
        ax.set_xlabel('Length (X)')
        ax.set_ylabel('Transverse (Y)')
        ax.set_zlabel('Vertical (Z)')
        ax.set_ylim(-10, 10)
        ax.set_zlim(-10, 10)
        
    plt.tight_layout()
    plt.show()

# --- MAIN EXECUTION (GLOBAL SCOPE) ---
# Because this code is NOT inside a function, variables will persist in Spyder.

if __name__ == "__main__":
    
    # 1. LOAD
    data = load_data()
    
    if data is not None:
        # 2. COMPUTE
        # Now returns S_matrix with ALL singular values
        freqs, S_mat, u_mats = compute_fdd(data, SAMPLING_FREQUENCY, N_PER_SEG)
        
        # Extract just the 1st singular value for peak picking
        s1_main = S_mat[:, 0]
        
        # 3. PICK PEAKS
        peaks = get_smart_peaks(freqs, s1_main)
        
        print(f"\n{'='*30}\nIDENTIFIED MODES\n{'='*30}")
        for i, p in enumerate(peaks):
            print(f"Mode {i+1}: {freqs[p]:.3f} Hz")
            
        # 4. PLOT SPECTRUM - UPDATED FOR MULTIPLE SVDs and NEW LIMITS
        plt.figure(figsize=(12, 6))
        
        # A) Plot the "Other" SVDs (3rd, 4th, ...) first so they are in background
        num_svds = S_mat.shape[1]
        if num_svds > 2:
            # Plot the first one with a label, the rest without label to avoid legend clutter
            plt.semilogy(freqs, S_mat[:, 2], 'k:', linewidth=0.8, alpha=0.4, label="Higher SVDs")
            for k in range(3, num_svds):
                plt.semilogy(freqs, S_mat[:, k], 'k:', linewidth=0.8, alpha=0.4)
        
        # B) Plot the 2nd SVD (Green, dashed)
        if num_svds > 1:
            plt.semilogy(freqs, S_mat[:, 1], 'g--', linewidth=1.5, alpha=0.8, label="2nd Singular Value")
            
        # C) Plot the 1st SVD (Blue, Solid - Main)
        plt.semilogy(freqs, s1_main, 'b-', linewidth=1.5, alpha=1.0, label="1st Singular Value")
        
        # Plot Peaks on top
        plt.plot(freqs[peaks], s1_main[peaks], 'rx', markersize=10, markeredgewidth=2, label="Identified Peaks")
        
        for i, p in enumerate(peaks):
            plt.annotate(f"{freqs[p]:.2f}", (freqs[p], s1_main[p]), xytext=(0, 10), 
                         textcoords='offset points', ha='center', color='red')
            
        plt.title("Frequency Domain Decomposition (Multiple SVDs)")
        plt.xlabel("Frequency (Hz)")
        plt.ylabel("Singular Values (dB)")
        
        # --- NEW LIMIT HERE ---
        # Adjusted to cover exactly the area of the first 6 modes
        plt.xlim(0, 13) 
        
        plt.grid(True, which='both', alpha=0.3)
        plt.legend(loc="upper right")
        plt.show()
        
        # 5. PLOT SHAPES
        plot_stick_model(u_mats, peaks, freqs)