%Example of Wave Distortion - wave clipping
%Steven Sahyun, University of Wisconsin - Whitewater, Nov. 12, 2021
clear;
clf;

bps = 16;       % bits per sample
sps = 16000;     % sample rate [samples/s]
freq = 300;       % frequency of low tone [Hz]
nsecs = 2;      % number of seconds of the audio file

% Set wave amplitude, noise, and nonlinearity values. 
% Amplitude values greater than clipAmp will be clipped.

amp = 2.5; % signal amplitude
clipAmp = 5; % signal clipping amplitude [max = 10]
noiseAmp = 2; % white noise amplitude [0 to + values]
nonlinear = 1% creates a power relationship to the wave amplitude [-0.49 to + values]

% Set Wave clipping values
waveMax = clipAmp;
waveMin = -clipAmp;

%Generate the wave
nsamples = sps*nsecs;
time = linspace(0, nsecs, nsamples);
wave = amp*sin(time*2*pi*freq)'; % a basic sine wave.

%Clipping Distortion
wave(wave>waveMax) = waveMax;
wave(wave<waveMin) = waveMin;


%noisy wave - add a white noise background
x = noiseAmp*randn(nsamples, 1);
fx = fft(x);
N = length(fx);
p = 0.5/2;
P = round(p*N);
fx(P+1:end)=0;
wavenoise = real(ifft(fx));
wave = wave+wavenoise;


%Nonlinear distortion
wave = clipAmp*real(wave.^(1+2*nonlinear))/clipAmp^(1+2*nonlinear);

%For showing the frequency spectrum of the wave
wavefft = fft(wave);

% create the amplitude window (creates the fade in and fade out when multiplied with the sine wave)
a = floor(.3*nsamples);
b = floor(.3*nsamples);
window = 1;
%window = [linspace(0,1,a), ones(1,nsamples-a-b), linspace(1,0,b)]';

%display the wavepattern and it's frequency spectrum
subplot(2, 1, 1);
plot(time, wave, 'LineWidth', 2);
axis([0, 0.01, -12-noiseAmp, 12+noiseAmp])
title('Sound Wave')
xlabel('Time (s)')

subplot(2,1,2);
xmax = 10*freq;
plot(8000*time, real(wavefft), 'LineWidth', 2)
ymax = max(real(wavefft));
ymin = min(real(wavefft));
axis([0, 2000, ymin, ymax]) % use if want to see change in frequency
%axis([0, xmax, ymin, ymax]) % use if want autoscale
title('Sound Wave Spectrum')
xlabel('Frequency (Hz)')
xticks(0:200:2000); %set to something easy to read for the frequency.
%xticks(0:freq:xmax); %autoscaled something easy to read for the frequency.

%Play audio of the sound wave
player = audioplayer(wave.*window, sps, bps);
play(player);