0% found this document useful (0 votes)
12 views32 pages

Python Signal Processing Techniques

Uploaded by

elisha immanuel
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
12 views32 pages

Python Signal Processing Techniques

Uploaded by

elisha immanuel
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd

EC23S01 NUMERICAL AND SIGNAL PROCESSING PRACTICE THROUGH

PYTHON LABORATORY
Getting Started with Python and using the Plot command

Continuous-Time Signal
(a) Plot a sine wave x(t)=sin⁡(2πft)x(t) = \sin(2\pi f t)x(t)=sin(2πft)

import numpy as np
import [Link] as plt
from scipy import signal
# Create a time array from 0 to 1 second, with 1000 samples
t = [Link](0, 1, 1000, endpoint=False)
# Create a sine wave signal
freq = 5 # 5 Hz
x = [Link](2 * [Link] * freq * t)
# Plot the signal
[Link](t, x)
[Link]("Sine Wave -5 Hz")
[Link]("Time [s]")
[Link]("Amplitude")
[Link](True)

[Link]()
(b) Square Wave

import numpy as np
import [Link] as plt
x_square = [Link](2 * [Link] * freq * t)
[Link](t, x_square)
[Link]("Square Wave -5 Hz")
[Link]("Time [s]")
[Link]("Amplitude")
[Link](True)
[Link]()

(c) Unit Step Function


import numpy as np
import [Link] as plt
from scipy import signal

x_square = [Link](2 * [Link] * freq * t)


t = [Link](-1, 1, 1000)
x = [Link](t, 1)
[Link](t, x)
[Link]("Unit Step Function")
1. Scaling a Signal
To multiply each sample of a signal x[n]x[n]x[n] by a constant aaa:
y[n]=a⋅x[n]y[n
2. Signal Addition
Adding two signals x[n] and y[n]
z[n]=x[n]+y[n]
3. Reversing a Signal
To reverse a signal x[n]→x[−n]
xrev[n]=x[−n] (represented by reversing the list)
4. Shifting a Signal
 Right Shift by k: y[n]=x[n−k]
 Left Shift by k: y[n]=x[n+k]

Procedure:
 Create a time vector n using range() or a list.
 Define discrete-time signals (step, ramp) using list comprehensions.
 Apply signal operations: scaling, addition, reversal, shifting.
 Use [Link]() to plot signals.
 Interpret the behavior of each signal and operation.

***************
import [Link] as plt
# 1. Creating basic signals using lists
n = list(range(0, 10)) # Time indices
# Unit step signal
unit_step = [1for i in n]
# Ramp signal
ramp = [i for i in n]
# 2. Scaling a signal
scaled_ramp = [2* x for x in ramp]
# 3. Adding two signals
combined_signal = [u + r for u, r in zip(unit_step, ramp)]
# 4. Reversing a signal
reversed_ramp = ramp[::-1]
# 5. Shifting a signal (Right shift by 2)
right_shifted = [0, 0] + ramp[:-2]
# Plotting all signals
[Link](figsize=(12, 8))
[Link](3, 2, 1)
[Link](n, unit_step)
[Link]("Unit Step Signal")
[Link]("n"); [Link]("Amplitude")
[Link](3, 2, 2)
[Link](n, ramp)
[Link]("Ramp Signal")
[Link](3, 2, 3)
[Link](n, scaled_ramp)
[Link]("Scaled Ramp Signal")
[Link](3, 2, 4)
[Link](n, combined_signal)
[Link]("Combined Signal (Step + Ramp)")
[Link](3, 2, 5)
[Link](n, reversed_ramp)
[Link]("Reversed Ramp Signal")
[Link](3, 2, 6)
[Link](n, right_shifted)
[Link]("Right Shifted Ramp Signal")
plt.tight_layout()
[Link]()
2. Load data from files and Plotting data

(a) sine wave import


(b) Square Wave

import numpy as np
import [Link] as plt
from scipy import signal
# Time vector
t = [Link](0, 1, 500, endpoint=False)
# Generate square wave
square_wave = [Link](2 * [Link] * 5 * t) # 5 Hz frequency
# Plotting
[Link](t, square_wave)
[Link]("5 Hz Square Wave")
[Link]("Time [s]")
[Link]("Amplitude")
[Link](True)
[Link]()
(c) Unit Step Function

import numpy as np
import [Link] as plt
# Time range
t = [Link](-5, 5, 1000)
# Unit step function
u = [Link](t >= 0, 1, 0)
# Plot
[Link](figsize=(8, 4))
[Link](t, u, label='Unit Step Function u(t)')
[Link]('Time')
[Link]('Amplitude')
[Link]('Unit Step Function')
[Link](True)
[Link]()
[Link]()
Exp. No.: 3 Getting Started with Date:
Lists
import math
import [Link] as plt
# Define n as a list
n = list(range(0, 10))
# Rectangular pulse x1[n] = 1 for n=0 to 4
x1 = list([1] * 5)
# Exponential decay signal x2[n] = exp(-0.5 * n)
x2 = [[Link](-0.5 * i) for i in n]
# Perform linear convolution manually using lists
m = len(x1)
n2 = len(x2)
conv_result = list([0] * (m + n2 - 1))
for i in range(m + n2 - 1):
for j in range(m):
if 0 <= i - j < n2:
conv_result[i] += x1[j] * x2[i - j]
# Plotting
[Link](figsize=(10, 6))
# x1[n]
[Link](3, 1, 1)
[Link](range(len(x1)), x1, linefmt='b-', markerfmt='bo', basefmt=" ", label='Rectangular Pulse')
[Link]()
# x2[n]
[Link](3, 1, 2)
[Link](range(len(x2)), x2, linefmt='g-', markerfmt='go', basefmt=" ", label='Exponential Decay')
[Link]()
# Convolution result
[Link](3, 1, 3)
[Link](range(len(conv_result)), conv_result, linefmt='r-', markerfmt='ro', basefmt="",
label='Convolution Result')
[Link]()
[Link]('n')
[Link](True)
plt.tight_layout()
[Link]()

***linear conv

import [Link] as plt


# Generate Input Signals using list()
x = list([1, 2, 3, 4])
h = list([1, 1, 1, -1])
# Lengths
m = len(x)
n = len(h)
# Zero Padding using list() and +
X = x + list([0]*n)
H = h + list([0]*m)
# Convolution Operation using lists
Y = list([0]*(m + n - 1))
for i in range(m + n - 1):
for j in range(m):
if 0 <= i - j < n + m:
Y[i] += X[j] * H[i - j]
# Plot results
[Link](figsize=(10, 6))
[Link](3, 1, 1)
[Link](range(len(x)), x, linefmt='b-', markerfmt='b^',
basefmt='k')
[Link]('n')
[Link]('x[n]')
[Link](True)
[Link](3, 1, 2)
[Link](range(len(h)), h, linefmt='m-', markerfmt='ms',
basefmt='k')
[Link]('n')
[Link]('h[n]')
[Link](True)
[Link](3, 1, 3)
[Link](range(len(Y)), Y, linefmt='r-', markerfmt='ro',
basefmt='k')
[Link]('n')
[Link]('Y[n]')
[Link](True)
plt.tight_layout()
[Link]()
circular
************
import numpy as np
import [Link] as plt
# Input sequences using list()
x = list([1, 2, 3, 4])
h = list([1, 1, 1, -1])
# Length for circular convolution = max(len(x), len(h))
N = max(len(x), len(h))
# Zero pad both sequences to length N using list() and +
x_padded = x + list([0] * (N - len(x)))
h_padded = h + list([0] * (N - len(h)))
# Perform FFT-based circular convolution
X = [Link](x_padded)
H = [Link](h_padded)
Y = [Link](X * H).real # Keep real part
# Convert result to list
Y = list(Y)
# Plot results
[Link](figsize=(10, 6))
[Link](3, 1, 1)
[Link](range(len(x)), x, linefmt='b-', markerfmt='b^',
basefmt='k')
[Link]('n')
[Link]('x[n]')
[Link](True)
[Link](3, 1, 2)
[Link](range(len(h)), h, linefmt='m-', markerfmt='ms',
basefmt='k')
[Link]('n')
[Link]('h[n]')
[Link](True)
[Link](3, 1, 3)
[Link](range(len(Y)), Y, linefmt='r-', markerfmt='ro',
basefmt='k')
[Link]('n')
[Link]('Y[n]')
[Link](True)
plt.tight_layout()
[Link]()
import numpy as np
import [Link] as plt
from [Link] import butter, filtfilt

# Sample input signal


x = [1, 2, 3, 4] # You can use any real or complex signal
N = len(x)
X = [Link](x) # Use numpy's FFT for efficiency and correctness

# Print the DFT result


print("DFT of x[n]:")
for i, val in enumerate(X):
print(f"X[{i}] = {val}")

# Parameters
fs_list = [10, 20, 50]
f, duration, noise_std = 5, 2, 0.5
cutoff = 3

# Define low-pass filter


def lowpass_filter(data, cutoff, fs):
b, a = butter(4, cutoff / (0.5 * fs), btype='low')
return filtfilt(b, a, data)

# High-resolution time axis


time_fine = [Link](0, duration, 1000)
clean_signal = [Link](2 * [Link] * f * time_fine)

# Plotting
[Link](figsize=(10, 8))
for i, fs in enumerate(fs_list):
# Coarse sampling for noisy signal
time = [Link](0, duration, int(fs * duration), endpoint=False)
noisy_signal = [Link](2 * [Link] * f * time) + [Link](0, noise_std, len(time))
filtered_signal = lowpass_filter(noisy_signal, cutoff, fs)

# Interpolate to high-resolution time base


noisy_interp = [Link](time_fine, time, noisy_signal)
filtered_interp = [Link](time_fine, time, filtered_signal)

# Plot
[Link](len(fs_list), 1, i + 1)
[Link](time_fine, noisy_interp, 'r--', alpha=0.5, label='Noisy')
[Link](time_fine, filtered_interp, 'b', label='Filtered')
[Link](time_fine, clean_signal, 'k', alpha=0.6, label='Clean')
[Link](f"Sampling Rate: {fs} Hz")
[Link]("Time (s)")
[Link]("Amplitude")
[Link](loc="upper right")

plt.tight_layout()
[Link]()

Exp. No.: 4 Getting started with for, Date:


If, While loops

Fourier Transform of Basic Signals


Write a Python program to compute the Discrete Fourier Transform (DFT) of a sinusoidal signal
x(t)=sin(2π10t). Use [Link] to perform the Fourier Transform and plot both the time-domain
signal and the magnitude spectrum. Discuss how the frequency content is reflected in the spectrum.

DFT of the signal using inbuilt function

import numpy as np
import [Link] as plt
t = [Link](0, 1,500)
x = [Link](2 * [Link] * 100 * t)+[Link](2 * [Link] * 200 * t)
X_f = [Link](x)
frequencies = [Link](len(t), d=t[1] - t[0])
#print(frequencies)
magnitude_spectrum = [Link](X_f)
#print(magnitude_spectrum)
[Link](figsize=(12, 6))
[Link](2, 1, 1)
[Link](t, x, label='x(t) = sin(2π10t)')
[Link]('Time (s)')
[Link]('Amplitude')
[Link]('Time-Domain Signal')
[Link](True)
[Link](2, 1, 2)
[Link](frequencies[:len(frequencies)//2],
magnitude_spectrum[:len(frequencies)//2])
[Link]('Frequency (Hz)')
[Link]('Magnitude')
[Link]('Magnitude Spectrum')
[Link](True)
plt.tight_layout()
[Link]()

Dft of the signal using formalue

x = [1, 2, 3, 4] # You can use any real or complex signal


N = len(x)
X = []
# DFT formula:
#X[k] = sum_n=0_to_N-1 x[n] * exp(-j*2*pi*k*n/N)
for k in range(N):
real = 0
imag = 0
for n in range(N):
angle = 2 * [Link] * k * n / N
real += x[n] * [Link](-angle)
imag += x[n] * [Link](-angle)
[Link](complex(real, imag))

#Print the DFT result


print("DFT of x[n]:")
for i, val in enumerate(X):
print(f"X[{i}] = {val}")

X[0] = (10+0j)
X[1] = (-2.0000000000000004+1.9999999999999996j)
X[2] = (-2-9.797174393178826e-16j)
X[3] = (-1.9999999999999982-2.000000000000001j)

Using for statement signal filtering

**
import numpy as np
import [Link] as plt
from [Link] import butter, filtfilt
# Parameters
fs_list = [10, 20, 50]
f, duration, noise_std = 5, 2, 0.5
cutoff = 3
# Define low-pass filter
def lowpass_filter(data, cutoff, fs):
b, a = butter(4, cutoff / (0.5 * fs), btype='low')
return filtfilt(b, a, data)
# High-resolution time axis
time_fine = [Link](0, duration, 1000)
clean_signal = [Link](2 * [Link] * f * time_fine)
# Plotting
[Link](figsize=(10, 8))
for i, fs in enumerate(fs_list):
# Coarse sampling for noisy signal
time = [Link](0, duration, int(fs * duration), endpoint=False)
noisy_signal = [Link](2 * [Link] * f * time) + [Link](0,
noise_std, len(time))
filtered_signal = lowpass_filter(noisy_signal, cutoff, fs)
# Interpolate to high-resolution time base
noisy_interp = [Link](time_fine, time, noisy_signal)
filtered_interp = [Link](time_fine, time, filtered_signal)
# Plot
[Link](len (fs_list), 1, i+1)
print(i)
print(fs_list)
[Link](time_fine, noisy_interp, 'r--', alpha=0.5, label='Noisy')
[Link](time_fine, filtered_interp, 'b', label='Filtered')
[Link](time_fine, clean_signal, 'k', alpha=0.6, label='Clean')
[Link](f"Sampling Rate: {fs} Hz")
[Link]("Time (s)")
[Link]("Amplitude")
[Link](loc="upper right")
plt.tight_layout()
[Link]()

# Original signal (can be any list of values)


signal = [1, 3, 7, 12, 5, -4, -10, 2]
# Clipping thresholds
min_val = -5
max_val = 10
# Clipped signal (using 'if')
clipped_signal = []
for sample in signal:
if sample > max_val:
clipped_signal.append(max_val)
elif sample < min_val:
clipped_signal.append(min_val)
else:
clipped_signal.append(sample)
print("Original Signal: ", signal)
print("Clipped Signal: ", clipped_signal)

while loop for signal normalizxation

signal = [4, 8, 15, 16, 23, 42]


# Find min and max
min_val = min(signal)
max_val = max(signal)
# Normalized signal (initialize)
normalized_signal = []
# Counter for while loop
i = 0
while i<len(signal):
# Apply normalization formula
norm = (signal[i] - min_val) / (max_val - min_val)
normalized_signal.append(norm)
i += 1
# Display results
print("Original Signal: ", signal)
print("Normalized Signal: ", normalized_signal)

Original Signal: [4, 8, 15, 16, 23, 42]


Normalized Signal: [0.0, 0.10526315789473684, 0.2894736842105263,
0.3157894736842105, 0.5, 1.0]

&&&&&&&&&&&&&&&&&&&&&&&&&
5/Getting started with files and arrays Aliasing in Real-World Signals

import numpy as np
import [Link] as plt
from [Link] import resample
# Original signal parameters
sr = 44100
t = [Link](0, 1, sr, endpoint=False)
original_signal = [Link](2 * [Link] * 440 * t) + 0.5 * [Link](2 * [Link] *
880 * t)
# Downsampling rates
downsampling_rates = [8000, 16000, 44100]
# Create figure
[Link](figsize=(12, 8))
for idx, rate in enumerate(downsampling_rates):
# Downsample the signal
downsampled_signal = resample(original_signal, rate)
# Compute FFT and frequency axis
frequencies = [Link](len(downsampled_signal), d=1/rate)
Y_f = [Link](downsampled_signal)
# Plot time-domain signal
[Link](3, 2, 2 * idx + 1)
[Link](downsampled_signal)
[Link](f'Downsampled Signal at {rate} Hz')
[Link]('Sample')
[Link]('Amplitude')
# Plot frequency-domain spectrum
[Link](3, 2, 2 * idx + 2)
[Link](frequencies[:len(frequencies)//2], [Link](Y_f)[:len(Y_f)//2])
[Link](f'Frequency Spectrum at {rate} Hz')
[Link]('Frequency (Hz)')
[Link]('Magnitude')
# Final layout adjustment
plt.tight_layout()
[Link]()
CODE: Using librosa
import librosa
import numpy as np
import [Link] as plt
from [Link] import Audio
# Load audio file (example trumpet sound)
signal, sr = [Link]([Link]('trumpet'))
# Function to plot frequency spectrum
def plot_spec(sig, sr, tle="Spectrum"):
D = [Link]([Link](sig))
freqs = [Link](len(D), 1/sr)

[Link](figsize=(8, 4))
[Link](freqs[:len(freqs)//2], D[:len(D)//2])
[Link](tle) # Fixed typo from 'tle'
[Link]("Frequency (Hz)")
[Link]("Amplitude")
[Link](True)
[Link]()
# Plot spectrum of original signal
plot_spec(signal, sr, "Original Signal Spectrum")
# Downsample rates
rates = [8000, 16000, 44100]
# Downsample, plot spectrum, and play audio
for r in rates:
ds_signal = [Link](y=signal, orig_sr=sr, target_sr=r)
plot_spec(ds_signal, r, f"Downsampled to {r} Hz")

import librosa
import numpy as np
import [Link] as plt
from [Link] import Audio
# Load audio file (example trumpet sound)
signal, sr = [Link]([Link]('trumpet'))
# Function to plot frequency spectrum
def plot_spec(sig, sr, tle="Spectrum"):
D = [Link]([Link](sig))
freqs = [Link](len(D), 1/sr)
[Link](figsize=(8, 4))
[Link](freqs[:len(freqs)//2], D[:len(D)//2])
[Link](tle) # Fixed typo from 'tle'
[Link]("Frequency (Hz)")
[Link]("Amplitude")
[Link](True)
[Link]()
# Plot spectrum of original signal
plot_spec(signal, sr, "Original Signal Spectrum")
# Downsample rates
rates = [8000, 16000, 44100]
# Downsample, plot spectrum, and play audio
for r in rates:
ds_signal = [Link](y=signal, orig_sr=sr, target_sr=r)
plot_spec(ds_signal, r, f"Downsampled to {r} Hz")
Power Spectral Density (PSD) - Using Numpy Array
Write a Python program to compute the Power Spectral Density (PSD) of a noisy signal consisting of
a 5 Hz sinusoid and Gaussian noise. Use the matplotlib library to plot the PSD and explain how the
noise and signal components appear in the frequency

import numpy as np
import [Link] as plt
from [Link] import welch
# Sampling settings
fs = 100 # Sampling frequency (Hz)
t_duration = 2 # Duration in seconds
f_signal = 5 # Frequency of the sine wave (Hz)
# Time axis - creates a 1D NumPy array of time points
t = [Link](0, t_duration, 1/fs)
# Signal + noise

# Creates a sine wave using the time array t


signal = [Link](2 * [Link] * f_signal * t)
# Generates Gaussian noise as a NumPy array of the same length as t
noise = [Link](0, 0.5, len(t))
# Element-wise addition of two arrays (signal + noise).
noisy_signal = signal + noise

# Power Spectral Density using Welch method


frequencies, psd = welch(noisy_signal, fs, nperseg=256)
# Plot
[Link](figsize=(10, 6))
# Plotting two arrays: frequency values vs. power values.
[Link](frequencies, psd, label='PSD of Noisy Signal')
[Link]('Power Spectral Density (PSD)')
[Link]('Frequency (Hz)')
[Link]('Power/Frequency (dB/Hz)')
[Link](True)
[Link]()
[Link](0, 30)
[Link](1e-5, 1e1)
plt.tight_layout()
[Link]()

import numpy as np
import [Link] as plt
from [Link] import welch
# Sampling settings
fs = 100 # Sampling frequency (Hz)
t_duration = 2 # Duration in seconds
f_signal = 5 # Frequency of the sine wave (Hz)
# Time axis - creates a 1D NumPy array of time points
t = [Link](0, t_duration, 1/fs)
# Signal + noise

# Creates a sine wave using the time array t


signal = [Link](2 * [Link] * f_signal * t)
# Generates Gaussian noise as a NumPy array of the same length as t
noise = [Link](0, 0.5, len(t))
# Element-wise addition of two arrays (signal + noise).
noisy_signal = signal + noise

# Power Spectral Density using Welch method


frequencies, psd = welch(noisy_signal, fs, nperseg=256)
# Plot
[Link](figsize=(10, 6))
# Plotting two arrays: frequency values vs. power values.
[Link](frequencies, psd, label='PSD of Noisy Signal')
[Link]('Power Spectral Density (PSD)')
[Link]('Frequency (Hz)')
[Link]('Power/Frequency (dB/Hz)')
[Link](True)
[Link]()
[Link](0, 30)
[Link](1e-5, 1e1)
plt.tight_layout()
[Link]()
Statistics using Python

import numpy as np
import [Link] as plt
from [Link] import skew, kurtosis

# Generate a noisy sine signal


fs = 1000
t = [Link](0, 1, fs, endpoint=False)
signal = [Link](2 * [Link] * 5 * t)
noise = [Link](0, 0.5, size=fs)
noisy_signal = signal + noise

# --- Statistics ---


mean_val = [Link](noisy_signal)
std_val = [Link](noisy_signal)
rms_val = [Link]([Link](noisy_signal**2))
skew_val = skew(noisy_signal)
kurt_val = kurtosis(noisy_signal)

# --- Output Statistics ---


print(f"Mean: {mean_val:.4f}")
print(f"Standard Deviation: {std_val:.4f}")
print(f"RMS: {rms_val:.4f}")
print(f"Skewness: {skew_val:.4f}")
print(f"Kurtosis: {kurt_val:.4f}")

# --- Plot Signal & Histogram ---


[Link](figsize=(12, 5))
[Link](1, 2, 1)
[Link](t, noisy_signal)
[Link]("Noisy Signal")
[Link]("Time (s)")
[Link]("Amplitude")

[Link](1, 2, 2)
[Link](noisy_signal, bins=50, density=True, color='skyblue', edgecolor='black')
[Link]("Signal Amplitude Histogram")
[Link]("Amplitude")
[Link]("Probability Density")

plt.tight_layout()
[Link]()
Linear combination of vectors and Computation of determinant, rank of a matrix

import numpy as np
# Define vectors
v1 = [Link]([1, 2])
v2 = [Link]([3, 4])
c1 = 2
c2 = -1

# Linear combination: v = c1*v1 + c2*v2


v = c1 * v1 + c2 * v2
print("Vector v1:", v1)
print("Vector v2:", v2)
print("Linear combination (2*v1 - 1*v2):", v)

#2[1 2 ] -1[3 4]
#2-3 4-4
#[-1 0]

import numpy as np
# Define matrix
A = [Link]([[2, 3, 1],
[4, 7, 7],
[6, 18, 22]])

#2(154-126)-3(88-42)+1(72-42)
#2(28) -138+30
#56-138+30
#-52

[1 2 3
142
265]

# Compute determinant
det = [Link](A)
# Compute rank
rank = [Link].matrix_rank(A)
print("Matrix A:\n", A)
print("Determinant of A:", round(det, 4))
print("Rank of A:", rank)
import numpy as np
def lu_decomposition(A):
n = len(A)
L = [Link]((n, n))
U = [Link]((n, n))
for i in range(n):
# Upper Triangular
for k in range(i, n):
sum_ = sum(L[i][j] * U[j][k] for j in range(i))
U[i][k] = A[i][k] - sum_
# Lower Triangular
for k in range(i, n):
if i == k:
L[i][i] = 1
else:
sum_ = sum(L[k][j] * U[j][i] for j in range(i))

L[k][i] = (A[k][i] - sum_) / U[i][i]


return L, U
def forward_substitution(L, b):
n = len(b)
y = [Link](n)
for i in range(n):
y[i] = b[i] - sum(L[i][j] * y[j] for j in range(i))
return y
def backward_substitution(U, y):
n = len(y)
x = [Link](n)
for i in reversed(range(n)):
x[i] = (y[i] - sum(U[i][j] * x[j] for j in range(i+1, n))) / U[i][i]
return x
# Main function
A = [Link]([[2, 3, 1],
[4, 7, 7],
[6, 18, 22]], dtype=float)
b = [Link]([1, 2, 3], dtype=float)
L, U = lu_decomposition(A)
y = forward_substitution(L, b)
x = backward_substitution(U, y)
print("Matrix L:\n", L)
print("Matrix U:\n", U)
print("Solution x:\n", x)

import numpy as np
def lu_decomposition(A):
n = len(A)
L = [Link]((n, n))
U = [Link]((n, n))
for i in range(n):
# Upper Triangular
for k in range(i, n):
sum_ = sum(L[i][j] * U[j][k] for j in range(i))
U[i][k] = A[i][k] - sum_
# Lower Triangular
for k in range(i, n):
if i == k:
L[i][i] = 1
else:
sum_ = sum(L[k][j] * U[j][i] for j in range(i))

L[k][i] = (A[k][i] - sum_) / U[i][i]


return L, U
def forward_substitution(L, b):
n = len(b)
y = [Link](n)
for i in range(n):
y[i] = b[i] - sum(L[i][j] * y[j] for j in range(i))
return y
def backward_substitution(U, y):
n = len(y)
x = [Link](n)
for i in reversed(range(n)):
x[i] = (y[i] - sum(U[i][j] * x[j] for j in range(i+1, n))) / U[i][i]
return x
# Main function
A = [Link]([[2, 3, 1], [4, 7, 7], [6, 18, 22]], dtype=float)
b = [Link]([1, 2, 3], dtype=float)
L, U = lu_decomposition(A)
y = forward_substitution(L, b)
x = backward_substitution(U, y)
print("Matrix L:\n", L)
print("Matrix U:\n", U)
print("Solution x:\n", x)

Generation of basic sequences using Python

import numpy as np
import [Link] as plt
n = [Link](-10, 11, 1)
# Unit Impulse
impulse = [Link](n == 0, 1, 0)
# Unit Step
step = [Link](n >= 0, 1, 0)
# Ramp
ramp = [Link](n >= 0, n, 0)
# Exponential
a = 0.8
expo = a ** n
# Sinusoidal
sin_seq = [Link](0.2 * [Link] * n)
# Plotting
[Link](figsize=(10, 10))
[Link](3, 2, 1)
[Link](n, impulse)
[Link]('Unit Impulse δ[n]')

[Link](3, 2, 2)
[Link](n, step)
[Link]('Unit Step u[n]')
[Link](3, 2, 3)
[Link](n, ramp)
[Link]('Ramp Sequence')
[Link](3, 2, 4)
[Link](n, expo)
[Link]('Exponential Sequence')
[Link](3, 2, 5)
[Link](n, sin_seq)
[Link]('Sinusoidal Sequence')
plt.tight_layout()
[Link]()

Generation of basic sequences using Python

import numpy as np
import [Link] as plt
n = [Link](-10, 11, 1)
# Unit Impulse
impulse = [Link](n == 0, 1, 0)
# Unit Step
step = [Link](n >= 0, 1, 0)
# Ramp
ramp = [Link](n >= 0, n, 0)
# Exponential
a = 0.8
expo = a ** n
# Sinusoidal
sin_seq = [Link](0.2 * [Link] * n)
# Plotting
[Link](figsize=(10, 10))
[Link](3, 2, 1)
[Link](n, impulse)
[Link]('Unit Impulse δ[n]')

[Link](3, 2, 2)
[Link](n, step)
[Link]('Unit Step u[n]')
[Link](3, 2, 3)
[Link](n, ramp)
[Link]('Ramp Sequence')
[Link](3, 2, 4)
[Link](n, expo)
[Link]('Exponential Sequence')
[Link](3, 2, 5)
[Link](n, sin_seq)
[Link]('Sinusoidal Sequence')
plt.tight_layout()
[Link]()
Spectral analysis of signals

import numpy as np
import [Link] as plt
# Generating a sample musical note signal
fs = 1100 # Sampling frequency (Hz)
duration = 2 # seconds
frequency = 440 # A4 note frequency (Hz)

t = [Link](0, duration, int(fs * duration), endpoint=False)


signal = [Link](2 * [Link] * frequency * t) + [Link](0, 1, len(t)) # Signal with noise
# Applying FFT
fft_result = [Link](signal)
freq = [Link]([Link][-1], d=1/fs)
# Plotting the spectrum
[Link](freq, [Link](fft_result))
[Link]('FFT of a Musical Note')
[Link]('Frequency (Hz)')
[Link]('Amplitude')
[Link]()

import numpy as np
import [Link] as plt
# Generating a sample musical note signal
fs = 1100 # Sampling frequency (Hz)
duration = 2 # seconds
frequency = 440 # A4 note frequency (Hz)

t = [Link](0, duration, int(fs * duration),


endpoint=False)
signal = [Link](2 * [Link] * frequency * t) +
[Link](0, 1, len(t)) # Signal with
noise
# Applying FFT
fft_result = [Link](signal)
freq = [Link]([Link][-1], d=1/fs)
# Plotting the spectrum
[Link](freq, [Link](fft_result))
[Link]('FFT of a Musical Note')
[Link]('Frequency (Hz)')
[Link]('Amplitude')
[Link]()

You might also like