1: Introduction to Digital Signal Processing

Learning Objectives

• Segue for DSP chapter.

Not only do we have analog signals --- signals that are real- or complex-valued functions of a continuous variable such as time or space --- we can define discrete/digital ones as well. Digital signals are sequences, functions defined only for the integers. We thus use the notation s(n) to denote a discrete-time one-dimensional signal such as a digital music recording and s(m,n) for a discrete-"time" two-dimensional signal like a photo taken with a digital camera. Sequences are fundamentally different than continuous-time signals. For example, continuity has no meaning for sequences.

Despite such fundamental differences, the theory underlying digital signal processing mirrors that of analog signals: Fourier transforms, linear filtering and linear systems parallel what previous chapters described. These similarities make it easy to understand the definitions and why we need them, but the similarities should not be construed as "analog wannabes." We will discover that digital signal processing is not an approximation to analog processing. We must explicitly worry about the fidelity of converting analog signals into digital ones. The music stored on CDs, the speech sent over digital cellular telephones, and the video carried by digital television all evidence that analog signals can be accurately converted to digital ones and back again.

The key reason why digital signal processing systems have a technological advantage today is the computer: computations, like the Fourier transform, can be performed quickly enough to be calculated as the signal is produced, and programmability means that the signal processing system can be easily changed. This flexibility has obvious appeal and has been widely accepted in the marketplace. Programmability means that we can perform signal processing operations impossible with analog systems (circuits). We will also discover that digital systems enjoy an algorithmic advantage that contributes to rapid processing speeds: Computations can be restructured in non-obvious ways to speed the processing. This flexibility comes at a price, a consequence of how computers work. How do computers perform signal processing?

Note

Taking a systems viewpoint for the moment, a system that produces its output as rapidly as the input arises is said to be a real-time system. All analog systems operate in real-time; digital ones that depend on a computer to perform system computations may or may not work in real-time. Clearly, we need real-time signal processing systems. Only recently have computers become fast enough to meet real-time requirements while performing non-trivial signal processing.

Digital signal processing has been applied to many different areas including medicine, telecommunication, military, and scientific to mention a few.  For example, the number of sun spots have been observed since the mid 1700's and analyzed to determine patterns of the frequency of their occurrence.   

Example

The following is an example of digital signal processing applied to scientific data analysis.  The sunspot data was obtained from the Solar Influences Data Analysis Center (SILSO) 

Figure \(\PageIndex{1}\): Sunspot Monthly Average Data.

Figure \(\PageIndex{2}\): Power Spectral Density of Sunspot Data.

The average number of sunspots that appeared on a monthly basis is displayed in the above graph, and the frequency analysis of these sunspots was obtained from power spectral density (PSD) analysis of this sunspot data as shown in the second graph.  Not surprisingly, there was a peak in the PSD analysis at 0.0938 cycles/year for a period of 10.67 years.  This corresponds to the observed cyclic variation of sunspots in the above graph.  The PSD analysis also reveals that there are two additional peaks: one at 0.0176 cycles/year with a period of 56.89 years and another at 0.1816 cycles/year for a 5.51 year period.  The ~57 year variation may be observed as the waxing and waning in the peak heights that occur every ~10.7 years corresponding to the completion of the sun's magnetic pole reversal.  The higher frequency peak that occurs at ~5.5 year cycle, is not noticeable in the sunspot data in the above figure which is consistent with its diminished peak magnitude in the PSD analysis relative to the other two peaks.

Python Code

 

'''
    Sunspot Number - Power Spectral Density
    Download Sunspot Data and Estimate Power Spectral Density
    Solar Influences Data Analysis Center (WDC-SILSO) [https://www.sidc.be/SILSO/home]
    Royal Observatory of Belgium, Brussels
    % Input data Monthly mean total sunspot number (1/1749 - 8/2024)
    Columns separated by one or more spaces
        Column 1-2: Gregorian calendar date
            - Year
            - Month
        Column 3: Date in fraction of year.
        Column 4: Monthly mean total sunspot number.
        Column 5: Monthly mean standard deviation of the input sunspot numbers.
        Column 6: Number of observations used to compute the monthly mean total
                sunspot number.
        Column 7: Definitive/provisional marker. A blank indicates that the value is
                definitive. A '*' indicates that the value is still provisional.
'''

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import sys

filename = "https://eng.libretexts.org/@api/deki/files/93479/SN_m_tot_V2.0.csv"
filename_url = "https://www.sidc.be/SILSO/DATA/SN_m_tot_V2.0.csv"   # Online source
'''
    Read Sunspot Data from local file or online source
    - File format: csv or txt
    - Columns separated by one or more spaces or semicolon
    - Use sep=';' for csv file or sep=r'\s+' for txt file with columns separated by one or more spaces
    - If the file is not found, it will try to read from the online source
    - If both fail, it will print an error message and set df to None
'''
# Attempt to read the local file first
try:
    df = pd.read_csv(filename,sep=';',index_col=False,header=None, \
            names=["Year","Month","Date","Sunspots","Sunspots_SD","Indicator","Provisional"])
except (FileNotFoundError, pd.errors.EmptyDataError):
    try:
        df = pd.read_csv(filename_url,sep=';',index_col=False,header=None, \
            names=["Year","Month","Date","Sunspots","Sunspots_SD","Indicator","Provisional"])
    except (FileNotFoundError, pd.errors.EmptyDataError):
        print("Error: Both file1.csv and file2.csv could not be loaded.")
        df = None # Or some default DataFrame
    else:
        print("Sunspot Data from Online source loaded successfully.\n")
except Exception as e:
    print(f"An unexpected error occurred: {e}")
else:
    print("Sunspot Data from local folder loaded successfully.\n")

if df is not None:
    # Proceed with operations on the DataFrame
    # Convert Year and Month to Date (YYYY-MM-DD format)
    df["Date"] = pd.to_datetime(df[["Year", "Month"]].assign(Day=1))
    print(df)   # Display DataFrame first and last 5 rows
else:
    print("No data to display. Please check the file paths or URLs.")
    sys.exit(1)

# Design the FIR lowpass filter & filter data
# Order of the filter
# Sampling frequency (cycles/year)
fs = 12     # Once per month (12 times per year)
# Cutoff frequency (cycles/year)
fc = 0.4    # Cutoff frequency cycles/year
order = 12
M1 = order + 1
# Lowpass Filter the Data
# Design the FIR lowpass filter using the window method
h1 = signal.firwin(M1, fc, window='hann', fs = fs)
y_filtered = signal.filtfilt(h1,[1],df["Sunspots"])

# Display Graph of Sunspots
plt.figure(figsize=(12,6))
plt.plot(df["Date"],df["Sunspots"], label='Monthly Averages')
plt.plot(df["Date"], y_filtered, 'r-', linewidth=1.5, label=f'Lowpass Filtered, fc = {fc} Cycles/Year')
plt.xlim(df["Date"].min(), df["Date"].max())
plt.grid()
plt.xlabel("Date (Year)")
plt.ylabel("Sunspot Number $S_n$")
plt.title("Numbers of Sunspots 1749 to Present (Monthly Averages), Lowpass Filtered with Raw Data")
plt.legend()

# Power Spectral Density Estimation
NFFT = 2048
f, Pxx = signal.welch(df["Sunspots"].to_numpy(),fs=fs,window='hann',nperseg=12*80, \
    noverlap=12*40,nfft=NFFT,detrend='linear',return_onesided=True,scaling='density')

plt.figure(figsize=(10,6))
plt.semilogy(f, Pxx)
plt.xlim(0,1)
plt.grid()
plt.xlabel('Frequency (Cycles/Year)')
plt.ylabel('Sunspots**2 - Year')
plt.title('Sunspot Power Spectral Density')
EndRange = int(0.25*NFFT/12)
peaks, _ = signal.find_peaks(Pxx[0:EndRange])
plt.text(0.2,5e4,rf'Peaks: {np.around(peaks*fs/NFFT, decimals=4)} Cycles/Year')
plt.text(0.2,2e4,rf'Periods: {np.around(np.reciprocal(peaks*fs/NFFT), decimals=2)} Years ')
plt.show()

Sunspot Data from Online source loaded successfully.

      Year  Month       Date  Sunspots  Sunspots_SD  Indicator  Provisional
0     1749      1 1749-01-01      96.7         -1.0         -1            1
1     1749      2 1749-02-01     104.3         -1.0         -1            1
2     1749      3 1749-03-01     116.7         -1.0         -1            1
3     1749      4 1749-04-01      92.8         -1.0         -1            1
4     1749      5 1749-05-01     141.7         -1.0         -1            1
...    ...    ...        ...       ...          ...        ...          ...
3312  2025      1 2025-01-01     137.0         23.3        670            0
3313  2025      2 2025-02-01     154.6         23.3        655            0
3314  2025      3 2025-03-01     134.2         20.4       1011            0
3315  2025      4 2025-04-01     140.6         17.6       1075            0
3316  2025      5 2025-05-01      79.2         10.6       1112            0

[3317 rows x 7 columns]