% Skript audio_1.m
clear;
% ------- Read audio file 
S = load('handel.mat');     % Data in S.y and sample freq in S.Fs
S

% ------- Write audiosignal in wav File
audiowrite('handel.wav',S.y,S.Fs);     % wav File
player = audioplayer(S.y,S.Fs);
play(player);
pause

figure(1);
plot((0:length(S.y)-1)/S.Fs, S.y);
La = axis;    axis([La(1:2),-1,1]);
title('handel');   xlabel('s');    grid on;

% ------- IIR-Filter
[B,A] = ellip(14,0.1,70,800/(S.Fs/2));   % IIR-Filter
nfft = 512;
HB = fft(B,nfft);
HA = fft(A,nfft);
Hiir = HB./HA;

figure(2);
plot((0:nfft-1)*S.Fs/nfft, 20*log10(abs(Hiir)),'k-','LineWidth',1);
title('Amplitude response of the IIR-Filter in dB')
xlabel('Hz');   grid on;

% -------- Signal plus noise
ny = length(S.y);
%noise = sqrt(std(S.y)^2/1000);       % Variance of noise
noise = sqrt(std(S.y)^2/10);       % Variance of noise

y1 = S.y + noise*randn(ny,1);         % Signal + noise

player = audioplayer(y1,S.Fs);
play(player);
pause;

figure(3);
plot((0:length(y1)-1)/S.Fs, y1,'k-','LineWidth',1);
hold on;
plot((0:length(S.y)-1)/S.Fs, S.y,'r-','LineWidth',1);

La = axis;    axis([La(1:2),-1,1]);
title('Signal and signal with noise');   xlabel('s');    grid on;

% ------ Filtering the Signal
y2 = filter(B,A,y1);

player = audioplayer(y2,S.Fs);
play(player);
pause

figure(4);
subplot(211), plot((0:length(y1)-1)/S.Fs, y1,'r-','LineWidth',1);
title('The signal with noise');

subplot(212), plot((0:length(y2)-1)/S.Fs, y2,'k-','LineWidth',1);
title('The filtered signal');

% ------- Power spectral density
[Py, w] = pwelch(S.y, hamming(512),[64],[512],S.Fs,'twosided');
[Py1, w] = pwelch(y1, hamming(512),[64],[512],S.Fs,'twosided');
[Py2, w] = pwelch(y2, hamming(512),[64],[512],S.Fs,'twosided');

figure(5);
subplot(311), plot(w, 10*log10(Py),'k-','LineWidth',1);
title('Poer spectral density of the original signal');
xlabel('Hz');  axis tight;

subplot(312), plot(w, 10*log10(Py1),'k-','LineWidth',1);
title('Poer spectral density of the signal with noise');
xlabel('Hz');  axis tight;

subplot(313), plot(w, 10*log10(Py2),'k-','LineWidth',1);
title('Poer spectral density of the filtered signal');
xlabel('Hz');  axis tight;




