classdef SignalLab % SignalLab Signal generation and analysis for ECE labs % % Provides methods for generating periodic waveforms, computing their % frequency spectra, and simulating measurement effects such as noise % and ADC quantization. Designed for use with the ADALM2000 (M2K). % % USAGE % s = SignalLab(amplitude, frequency, samplingRate) % s = SignalLab(amplitude, frequency, samplingRate, phase) % % CONFIGURATION % s = s.setWave("square") % change waveform type (must re-assign) % s = s.setWave("square", 75) % square wave with 75% duty cycle % s.info() % display all current settings % % GENERATION % [t, x] = s.generate(duration) % [t, x] = s.generateHarmonics(N, duration) % [t, x] = s.generateNoisy(duration, noiseAmp) % % ANALYSIS % A = s.measureAmplitude(x) % f = s.measureFrequency(x) % snr_db = s.calculateSNR(x) % xf = s.filterSignal(x, fLow, fHigh) % xq = s.quantizeSignal(x, range, bits) % % SPECTRAL ANALYSIS % [fDom, spectrum, freqAxis] = s.computeSpectrum(x) % % PLOTTING % s.plotSignal(t, x) % s.plotSpectrum(x) % s.visualizeSNR(x) % s.plotQuantization(x, range, bits) % % IMPORTANT: SignalLab is a value class. Methods that modify settings % must be re-assigned: s = s.setWave("square") % % Version 1.0 | ECE Emerge Lab properties amplitude % Peak signal amplitude (V) frequency % Fundamental frequency (Hz) phase % Phase offset (radians), default 0 samplingRate % Sampling rate (Hz) waveType % "sine" | "square" | "triangle" | "sawtooth" dutyCycle % Square wave duty cycle (0–100%), default 50 end % Maximum amplitude before a warning is issued (M2K input range limit) properties (Constant, Access = private) MAX_AMPLITUDE = 25 % V — M2K maximum input range MAX_FREQUENCY = 1e6 % Hz — practical upper limit for lab use MAX_SAMPLING_RATE = 1e8 % Hz — practical upper limit for lab use end methods % ----------------------------------------------------------------- % Constructor % ----------------------------------------------------------------- function obj = SignalLab(amplitude, frequency, samplingRate, phase) % SignalLab Create a signal generator/analyser object % % s = SignalLab(amplitude, frequency, samplingRate) % s = SignalLab(amplitude, frequency, samplingRate, phase) % % amplitude - Peak amplitude in volts % frequency - Frequency in Hz (must be > 0) % samplingRate - Sampling rate in Hz (must be > 0) % phase - Phase offset in radians (default: 0) % % Warnings are issued for: % - Zero amplitude (signal will be all zeros) % - |amplitude| > 25 V (exceeds M2K input range) % - |phase| > 2*pi (outside standard range) % - samplingRate < 2*frequency (violates Nyquist criterion) narginchk(3, 4); if nargin < 4 phase = 0; end % --- Type and basic range checks --- if ~isnumeric(amplitude) || ~isscalar(amplitude) error("SignalLab:invalidInput", "Amplitude must be a numeric scalar."); end if ~isnumeric(frequency) || ~isscalar(frequency) || frequency <= 0 error("SignalLab:invalidInput", "Frequency must be a positive numeric scalar."); end if frequency > SignalLab.MAX_FREQUENCY error("SignalLab:invalidInput", ... "Frequency (%.4g Hz) exceeds the practical limit of %.4g Hz.", ... frequency, SignalLab.MAX_FREQUENCY); end if ~isnumeric(samplingRate) || ~isscalar(samplingRate) || samplingRate <= 0 error("SignalLab:invalidInput", "Sampling rate must be a positive numeric scalar."); end if samplingRate > SignalLab.MAX_SAMPLING_RATE error("SignalLab:invalidInput", ... "Sampling rate (%.4g Hz) exceeds the practical limit of %.4g Hz.", ... samplingRate, SignalLab.MAX_SAMPLING_RATE); end if ~isnumeric(phase) || ~isscalar(phase) error("SignalLab:invalidInput", "Phase must be a numeric scalar."); end % --- Advisory warnings --- if amplitude == 0 warning("SignalLab:zeroAmplitude", ... "Amplitude is zero — the generated signal will be all zeros."); end if abs(amplitude) > SignalLab.MAX_AMPLITUDE warning("SignalLab:amplitudeRange", ... "Amplitude (%.4g V) exceeds the M2K input range of ±%.4g V.", ... amplitude, SignalLab.MAX_AMPLITUDE); end if abs(phase) > 2*pi warning("SignalLab:phaseRange", ... "Phase (%.4g rad) is outside the standard range of ±2*pi. " + ... "Consider using a value in [-2*pi, 2*pi].", phase); end if samplingRate < 2 * frequency warning("SignalLab:nyquist", ... "Sampling rate (%.4g Hz) is below the Nyquist rate (%.4g Hz). " + ... "The signal will be aliased.", ... samplingRate, 2*frequency); end obj.amplitude = amplitude; obj.frequency = frequency; obj.samplingRate = samplingRate; obj.phase = phase; obj.waveType = "sine"; obj.dutyCycle = 50; end % ----------------------------------------------------------------- % Configuration % ----------------------------------------------------------------- function obj = setWave(obj, waveType, dutyCycle) % setWave Change the waveform type % % s = s.setWave("sine") % s = s.setWave("square") % s = s.setWave("square", 75) % 75% duty cycle % % Supported types: "sine", "square", "triangle", "sawtooth" % % Note: duty cycles of 0% and 100% produce a constant DC % signal. A warning is issued if these values are used. validTypes = ["sine", "square", "triangle", "sawtooth"]; waveType = lower(string(waveType)); % accept char or string input if ~any(waveType == validTypes) error("SignalLab:invalidWaveType", ... "Wave type must be: ""sine"", ""square"", ""triangle"", or ""sawtooth""."); end obj.waveType = waveType; if nargin >= 3 if waveType ~= "square" warning("SignalLab:dutyCycleIgnored", ... "Duty cycle only applies to square waves — ignored for ""%s"".", ... waveType); else if ~isnumeric(dutyCycle) || ~isscalar(dutyCycle) || ... dutyCycle < 0 || dutyCycle > 100 error("SignalLab:invalidInput", ... "Duty cycle must be a numeric scalar between 0 and 100."); end if dutyCycle == 0 || dutyCycle == 100 warning("SignalLab:dutyCycleConstant", ... "Duty cycle of %g%% produces a constant DC signal.", ... dutyCycle); end obj.dutyCycle = dutyCycle; end end end function info(obj) % info Display all current signal properties % % s.info() % % Also shows the Nyquist frequency and the frequency % resolution achievable for a 1-second signal. nyquist = obj.samplingRate / 2; fprintf("--- SignalLab ---\n"); fprintf(" Amplitude : %g V\n", obj.amplitude); fprintf(" Frequency : %g Hz\n", obj.frequency); fprintf(" Phase : %g rad\n", obj.phase); fprintf(" Sampling Rate: %g Hz\n", obj.samplingRate); fprintf(" Nyquist Freq : %g Hz\n", nyquist); fprintf(" Wave Type : %s\n", obj.waveType); if obj.waveType == "square" fprintf(" Duty Cycle : %g%%\n", obj.dutyCycle); end if obj.samplingRate < 2 * obj.frequency fprintf(" *** WARNING: sampling rate is below the Nyquist rate ***\n"); end end % ----------------------------------------------------------------- % Signal generation % ----------------------------------------------------------------- function [t, x] = generate(obj, duration) % generate Generate the current waveform % % [t, x] = s.generate(duration) % % duration - Signal duration in seconds (positive scalar) % t - Time vector (s) % x - Signal samples (V) % % The duration is automatically adjusted to the nearest whole % number of complete periods to avoid spectral leakage. % % For sine waves: an error is raised if samplingRate < 2*frequency. % For non-sinusoidal waveforms: same error applies; additionally % a warning is issued if samplingRate < 10*frequency, because % ideal square/triangle/sawtooth waves have infinite bandwidth % and harmonics above Nyquist will fold back into the spectrum. obj.validateDuration(duration); obj.checkNyquist(); t = obj.timeVector(duration); tau = mod(2*pi*obj.frequency*t + obj.phase, 2*pi) / (2*pi); switch obj.waveType case "sine" x = obj.amplitude * sin(2*pi*obj.frequency*t + obj.phase); case "square" % +amplitude for first dutyCycle%, -amplitude for remainder x = obj.amplitude * (2*(tau < obj.dutyCycle/100) - 1); case "triangle" % Rises from -A to +A in first half period, back in second x = obj.amplitude * (1 - 4*abs(tau - 0.5)); case "sawtooth" % Rises linearly from -A to +A over one full period x = obj.amplitude * (2*tau - 1); end end function [t, x] = generateHarmonics(obj, N, duration) % generateHarmonics Fourier series approximation up to the Nth harmonic % % [t, x] = s.generateHarmonics(N, duration) % % N - Highest harmonic number to include (positive integer) % square/triangle : odd harmonics up to N (1, 3, 5, …, ≤N) % sawtooth : all harmonics up to N (1, 2, 3, …, N) % sine : N is ignored (sine has no harmonics) % duration - Signal duration in seconds % % Fourier series coefficients: % Square : x(t) = (4A/pi) * sum_{odd n ≤ N} sin(2*pi*n*f*t) / n % Triangle: x(t) = (8A/pi^2) * sum_{odd n ≤ N} (-1)^k * sin(2*pi*n*f*t) / n^2 % Sawtooth: x(t) = (2A/pi) * sum_{n=1}^{N} (-1)^(n+1) * sin(2*pi*n*f*t) / n % % A warning is issued if any harmonic frequency exceeds the % Nyquist limit — those harmonics would be aliased. if ~isnumeric(N) || ~isscalar(N) || N < 1 || N ~= floor(N) error("SignalLab:invalidInput", "N must be a positive integer."); end obj.validateDuration(duration); obj.checkNyquist(); % Check that highest harmonic does not exceed Nyquist highestFreq = N * obj.frequency; if highestFreq > obj.samplingRate / 2 warning("SignalLab:harmonicAliasing", ... "Harmonic %d (%.4g Hz) exceeds the Nyquist frequency (%.4g Hz) " + ... "and will be aliased. Reduce N or increase the sampling rate.", ... N, highestFreq, obj.samplingRate/2); end if obj.waveType == "sine" && N > 1 warning("SignalLab:harmonicsIgnored", ... "A sine wave has no harmonics above the fundamental. " + ... "N = %d is ignored — only the fundamental is generated.", N); end t = obj.timeVector(duration); x = zeros(size(t)); switch obj.waveType case "sine" x = obj.amplitude * sin(2*pi*obj.frequency*t + obj.phase); case "square" for n = 1 : 2 : N % odd harmonics: 1, 3, 5, …, ≤N x = x + sin(2*pi*n*obj.frequency*t + obj.phase) / n; end x = (4*obj.amplitude / pi) * x; case "triangle" for n = 1 : 2 : N % odd harmonics: 1, 3, 5, …, ≤N k = (n - 1) / 2; x = x + (-1)^k * sin(2*pi*n*obj.frequency*t + obj.phase) / n^2; end x = (8*obj.amplitude / pi^2) * x; case "sawtooth" for n = 1 : N % all harmonics: 1, 2, 3, …, N x = x + (-1)^(n+1) * sin(2*pi*n*obj.frequency*t + obj.phase) / n; end x = (2*obj.amplitude / pi) * x; end end function [t, x] = generateNoisy(obj, duration, noiseAmp) % generateNoisy Generate current waveform with added Gaussian noise % % [t, x] = s.generateNoisy(duration) % [t, x] = s.generateNoisy(duration, noiseAmp) % % noiseAmp - Standard deviation of the noise in volts (default: 0.1) % Must be a non-negative scalar. if nargin < 3 noiseAmp = 0.1; end if ~isnumeric(noiseAmp) || ~isscalar(noiseAmp) || noiseAmp < 0 error("SignalLab:invalidInput", ... "Noise amplitude must be a non-negative numeric scalar."); end [t, x] = obj.generate(duration); x = x + noiseAmp * randn(size(x)); end % ----------------------------------------------------------------- % Signal analysis % ----------------------------------------------------------------- function A = measureAmplitude(~, x) % measureAmplitude Estimate amplitude as half the peak-to-peak value % % A = s.measureAmplitude(x) % % x - Signal vector (numeric, at least 2 samples) % % Returns the peak amplitude in volts. SignalLab.validateSignalInput(x, "measureAmplitude"); A = (max(x) - min(x)) / 2; end function f = measureFrequency(obj, x) % measureFrequency Estimate the dominant frequency via FFT % % f = s.measureFrequency(x) % % x - Signal vector (numeric, at least 4 samples) % % Returns the frequency in Hz of the largest spectral peak, % excluding DC. SignalLab.validateSignalInput(x, "measureFrequency"); [f, ~, ~] = obj.computeSpectrum(x); end function snr_db = calculateSNR(obj, x) % calculateSNR Estimate signal-to-noise ratio in dB % % snr_db = s.calculateSNR(x) % % x - Signal vector (numeric, at least 4 samples) % % Signal power : spectral energy within ±5% of the dominant % frequency (minimum bandwidth 5 Hz) % Noise power : all other spectral energy, excluding DC % % Returns Inf if no noise is present. SignalLab.validateSignalInput(x, "calculateSNR"); [fDom, spectrum, freqAxis] = obj.computeSpectrum(x); bw = max(5, 0.05 * fDom); inBand = abs(freqAxis - fDom) <= bw; outBand = ~inBand & (freqAxis > 0); signalPower = sum(spectrum(inBand).^2); noisePower = sum(spectrum(outBand).^2); if noisePower == 0 snr_db = Inf; else snr_db = 10 * log10(signalPower / noisePower); end end function xf = filterSignal(obj, x, fLow, fHigh) % filterSignal Ideal rectangular bandpass filter % % xf = s.filterSignal(x, fLow, fHigh) % % x - Signal vector (numeric) % fLow, fHigh - Lower and upper cutoff frequencies in Hz % Must satisfy: 0 ≤ fLow < fHigh % % Frequencies outside [fLow, fHigh] are zeroed in the FFT % domain. Returns the filtered signal at the same length as x. % % A warning is issued if fHigh exceeds the Nyquist frequency % since no signal energy exists above Nyquist. SignalLab.validateSignalInput(x, "filterSignal"); if ~isnumeric(fLow) || ~isscalar(fLow) error("SignalLab:invalidInput", "fLow must be a numeric scalar."); end if ~isnumeric(fHigh) || ~isscalar(fHigh) error("SignalLab:invalidInput", "fHigh must be a numeric scalar."); end if fLow < 0 error("SignalLab:invalidInput", "fLow must be non-negative."); end if fLow >= fHigh error("SignalLab:invalidInput", "fLow must be less than fHigh."); end if fHigh > obj.samplingRate / 2 warning("SignalLab:nyquist", ... "fHigh (%.4g Hz) exceeds the Nyquist frequency (%.4g Hz). " + ... "No signal energy exists above Nyquist.", ... fHigh, obj.samplingRate/2); end x = x(:); N = length(x); Y = fft(x); freq = obj.samplingRate * (0:N-1) / N; % Build bandpass mask; mirror for the negative-frequency half mask = zeros(N, 1); for i = 1:N f = freq(i); fNeg = obj.samplingRate - f; if (f >= fLow && f <= fHigh) || (fNeg >= fLow && fNeg <= fHigh) mask(i) = 1; end end xf = real(ifft(Y .* mask)); end function xq = quantizeSignal(obj, x, range, bits) % quantizeSignal Simulate ADC quantization % % xq = s.quantizeSignal(x) % xq = s.quantizeSignal(x, range) % xq = s.quantizeSignal(x, range, bits) % % x - Input signal vector (V) % range - Full-scale ADC input range in volts % Default: 5 V (±2.5 V) — M2K low-voltage range % Use 50 V for the M2K high-voltage range (±25 V) % bits - ADC resolution in bits (default: 12 — M2K) % % Samples outside ±range/2 are clipped before quantization. % LSB size = range / 2^bits % % M2K ADC input ranges: % ±2.5 V (range = 5 V) — use for signals up to ±2.5 V % ±25 V (range = 50 V) — use for signals up to ±25 V SignalLab.validateSignalInput(x, "quantizeSignal"); if nargin < 3, range = 5; end if nargin < 4, bits = 12; end if ~isnumeric(range) || ~isscalar(range) || range <= 0 error("SignalLab:invalidInput", ... "Range must be a positive numeric scalar (full-scale volts)."); end if ~isnumeric(bits) || ~isscalar(bits) || bits < 1 || bits ~= floor(bits) error("SignalLab:invalidInput", "Bits must be a positive integer."); end nLevels = 2^bits; lsb = range / nLevels; % Clip to ADC input range xq = max(-range/2, min(range/2, x)); % Round to nearest quantization level xq = round(xq / lsb) * lsb; % Clamp upper edge to prevent rounding overshoot beyond +range/2 xq = min(xq, range/2 - lsb); end % ----------------------------------------------------------------- % Plotting % ----------------------------------------------------------------- function plotSignal(~, t, x) % plotSignal Plot the signal in the time domain % % s.plotSignal(t, x) % % t - Time vector (s) % x - Signal vector (V); must be the same length as t if ~isnumeric(t) || ~isvector(t) error("SignalLab:invalidInput", "t must be a numeric vector."); end if ~isnumeric(x) || ~isvector(x) error("SignalLab:invalidInput", "x must be a numeric vector."); end if length(t) ~= length(x) error("SignalLab:invalidInput", ... "t and x must be the same length (t: %d, x: %d).", ... length(t), length(x)); end plot(t, x, "b-", "LineWidth", 1.5); xlabel("Time (s)", "FontSize", 12); ylabel("Amplitude (V)", "FontSize", 12); title("Time Domain Signal", "FontSize", 14); grid on; box on; end function plotSpectrum(obj, x) % plotSpectrum Plot the single-sided magnitude spectrum in dB % % s.plotSpectrum(x) % % x - Signal vector (numeric, at least 4 samples) % % The FFT is computed internally — no prior setup needed. % The spectrum is scaled so that peak values equal the true % signal amplitudes before conversion to dB. SignalLab.validateSignalInput(x, "plotSpectrum"); [~, spectrum, freqAxis] = obj.computeSpectrum(x); plot(freqAxis, 20*log10(spectrum + eps), "r-", "LineWidth", 1.5); xlabel("Frequency (Hz)", "FontSize", 12); ylabel("Magnitude (dB)", "FontSize", 12); title("Frequency Spectrum", "FontSize", 14); grid on; box on; end function visualizeSNR(obj, x) % visualizeSNR Two-panel figure: time domain + spectrum with signal band % % s.visualizeSNR(x) % % x - Signal vector (numeric, at least 4 samples) % % The signal band (used for SNR calculation) is highlighted % in red on the spectrum panel. SignalLab.validateSignalInput(x, "visualizeSNR"); [fDom, spectrum, freqAxis] = obj.computeSpectrum(x); snr_db = obj.calculateSNR(x); bw = max(5, 0.05 * fDom); inBand = abs(freqAxis - fDom) <= bw; t = (0 : length(x)-1) / obj.samplingRate; figure("Position", [100, 100, 800, 600]); subplot(2, 1, 1); plot(t, x, "b-", "LineWidth", 1.2); xlabel("Time (s)"); ylabel("Amplitude (V)"); title("Time Domain Signal"); grid on; subplot(2, 1, 2); plot(freqAxis, spectrum, "b-", "LineWidth", 1.2); hold on; plot(freqAxis(inBand), spectrum(inBand), "r-", "LineWidth", 2.5); hold off; xlabel("Frequency (Hz)"); ylabel("Magnitude"); title(sprintf("Frequency Spectrum — SNR = %.1f dB", snr_db)); legend("Full Spectrum", "Signal Band", "Location", "northeast"); grid on; end function plotQuantization(obj, x, range, bits) % plotQuantization Visualise the effect of ADC quantization % % s.plotQuantization(x) % s.plotQuantization(x, range) % s.plotQuantization(x, range, bits) % % x - Signal vector (numeric, at least 4 samples) % range - Full-scale ADC range in volts (default: 5 V = ±2.5 V) % bits - ADC resolution in bits (default: 12) % % Three-panel figure: % Top — original vs quantized; clipped samples marked in red; % ADC limits shown as dotted lines % Middle — quantization error with ±LSB/2 bounds % Bottom — magnitude spectrum of original vs quantized; % theoretical noise floor shown in title SignalLab.validateSignalInput(x, "plotQuantization"); if nargin < 3, range = 5; end if nargin < 4, bits = 12; end xq = obj.quantizeSignal(x, range, bits); err = xq - x; t = (0 : length(x)-1) / obj.samplingRate; lsb = range / 2^bits; clipped = abs(x) > range/2; figure("Position", [100, 100, 900, 750]); % --- Top: original vs quantized, ADC limits and clipping --- subplot(3, 1, 1); plot(t, x, "b-", "LineWidth", 1.2); hold on; plot(t, xq, "r--", "LineWidth", 1.0); if any(clipped) plot(t(clipped), x(clipped), "ro", "MarkerSize", 4, ... "DisplayName", "Clipped"); end yline( range/2, "k:", "LineWidth", 1.0); yline(-range/2, "k:", "LineWidth", 1.0); hold off; xlabel("Time (s)"); ylabel("Amplitude (V)"); title(sprintf( ... "Original vs Quantized | %d-bit | range = \\pm%.4g V | LSB = %.4g mV", ... bits, range/2, lsb*1000)); if any(clipped) legend("Original", "Quantized", "Clipped", "Location", "northeast"); else legend("Original", "Quantized", "Location", "northeast"); end grid on; % --- Middle: quantization error --- subplot(3, 1, 2); plot(t, err * 1000, "k-", "LineWidth", 1.0); yline( lsb/2 * 1000, "r--", "+LSB/2", "LineWidth", 1.0); yline(-lsb/2 * 1000, "r--", "-LSB/2", "LineWidth", 1.0); xlabel("Time (s)"); ylabel("Error (mV)"); title(sprintf("Quantization Error | bounded by \\pm%.4g mV", lsb/2*1000)); grid on; % --- Bottom: spectrum comparison --- subplot(3, 1, 3); [~, spec_orig, freq] = obj.computeSpectrum(x); [~, spec_q, ~ ] = obj.computeSpectrum(xq); plot(freq, 20*log10(spec_orig + eps), "b-", "LineWidth", 1.2); hold on; plot(freq, 20*log10(spec_q + eps), "r--", "LineWidth", 1.2); hold off; xlabel("Frequency (Hz)"); ylabel("Magnitude (dB)"); title(sprintf( ... "Spectrum: Original vs Quantized | theoretical noise floor \\approx %.1f dB", ... -(6.02*bits + 1.76))); legend("Original", "Quantized", "Location", "northeast"); grid on; end end % ===================================================================== % Spectral analysis % ===================================================================== methods function [fDom, spectrum, freqAxis] = computeSpectrum(obj, x) % computeSpectrum Single-sided magnitude spectrum and dominant frequency % % [fDom, spectrum, freqAxis] = s.computeSpectrum(x) % % x - Signal vector (numeric, at least 4 samples) % fDom - Dominant frequency (Hz): frequency of the largest % spectral peak, excluding DC % spectrum - Single-sided magnitude spectrum (V); DC is removed % before the FFT and amplitudes are scaled so peaks % equal true signal amplitudes (e.g. a 1 V sine gives % a peak of 1 V at f0) % freqAxis - Frequency axis (Hz), length floor(N/2)+1, % with bin spacing Fs/N % % Example — extract the amplitude at the 3rd harmonic: % [~, spectrum, freqAxis] = s.computeSpectrum(x); % [~, idx] = min(abs(freqAxis - 3*s.frequency)); % amp3 = spectrum(idx); SignalLab.validateSignalInput(x, "computeSpectrum"); x = x(:) - mean(x); N = length(x); if N < 4 error("SignalLab:tooShort", ... "Signal must contain at least 4 samples for spectral analysis."); end Xmag = abs(fft(x) / N); nUniq = floor(N/2) + 1; spectrum = Xmag(1:nUniq); spectrum(2:end-1) = 2 * spectrum(2:end-1); freqAxis = (0 : nUniq-1) * (obj.samplingRate / N); % Dominant frequency: largest peak excluding DC (index 1) [~, idx] = max(spectrum(2:end)); fDom = freqAxis(idx + 1); end end % ===================================================================== % Private helpers % ===================================================================== methods (Access = private) function t = timeVector(obj, duration) % Snap to the nearest whole number of complete periods of the % fundamental frequency. This prevents spectral leakage in FFT % analysis without requiring the student to choose duration % carefully. numPeriods = max(1, round(duration * obj.frequency)); N = round(numPeriods * obj.samplingRate / obj.frequency); t = (0 : N-1) / obj.samplingRate; end function validateDuration(~, duration) % Validate that duration is a positive numeric scalar. if ~isnumeric(duration) || ~isscalar(duration) || duration <= 0 error("SignalLab:invalidInput", ... "Duration must be a positive numeric scalar (seconds)."); end end function checkNyquist(obj) % Raise an error or warning based on the Nyquist criterion. % % For ALL waveforms: % Error if samplingRate < 2 * frequency. % The fundamental cannot be represented at all. % % For non-sinusoidal waveforms (square, triangle, sawtooth): % Warning if samplingRate < 10 * frequency. % Ideal non-sinusoidal waveforms have infinite bandwidth and % are never truly bandlimited. The threshold of 10× ensures % at least the 1st through 4th odd harmonics (1, 3, 5, 7 ×f0) % fall below Nyquist. Below this threshold the waveform shape % will be visibly distorted by harmonic aliasing. if obj.samplingRate < 2 * obj.frequency error("SignalLab:nyquist", ... "Cannot generate signal: sampling rate (%.4g Hz) is below " + ... "the Nyquist rate (%.4g Hz) for frequency %.4g Hz. " + ... "Increase samplingRate or decrease frequency.", ... obj.samplingRate, 2*obj.frequency, obj.frequency); end if obj.waveType ~= "sine" && obj.samplingRate < 10 * obj.frequency highestHarmonic = floor(obj.samplingRate / (2 * obj.frequency)); warning("SignalLab:harmonicAliasing", ... "%s wave: sampling rate (%.4g Hz) is less than 10× the " + ... "fundamental (%.4g Hz). Only harmonics up to %d are below " + ... "the Nyquist frequency (%.4g Hz); higher harmonics will " + ... "fold back into the spectrum and distort the waveform shape. " + ... "Consider increasing the sampling rate.", ... obj.waveType, obj.samplingRate, obj.frequency, ... highestHarmonic, obj.samplingRate / 2); end end end % ===================================================================== % Static helpers % ===================================================================== methods (Static, Access = private) function validateSignalInput(x, callerName) % Validate that x is a non-empty numeric vector with at least % 2 samples. Raises a descriptive error referencing the caller. if ~isnumeric(x) || ~isvector(x) error("SignalLab:invalidInput", ... "%s: x must be a numeric vector.", callerName); end if length(x) < 2 error("SignalLab:invalidInput", ... "%s: signal must contain at least 2 samples.", callerName); end end end end