# Need help to find true Signal to noise ratio (SNR) for acceleration signal

8 views (last 30 days)
hanna Eid on 22 May 2022
Answered: Abhimenyu on 8 Dec 2023
Question:
I am quite new to matlab, and I have gotten the task to find the true Signal to noise ratio for the original acceleration signal. For the signal to noise ratio (SNR1) for the noisy acceleration signal, we got a value of 6.93.
So what I am struggling with is how to find the sum harmonics I need to use the code I was thinking. Does anybody know how to do that? or is it another better way?
file_path=uigetdir;
cd(file_path);
files=dir("acc_data*.mat");
Fig=1
Fs=100; %sampling frequency
F= Fs/2;
for i=1:length(files)
N=size(acc,2);
Win=10; %why 10 as a window size?
win_center=round(Win/2);
PSD_acc=[];Freq_acc=[];PSD_acc1=[];Freq_acc1=[];
for n=1:N
sd_acc=std(acc(indx,n)); %Standard deviation of acc indx
m_acc=mean(acc(indx,:)); %mean ...
acc_noise= (0.3*sd_acc)*randn(size(indx,1),n); %task 4:creating noisy version of signal
NOISYsignal=acc_noise+acc(indx,n); %Add noise with amplitude 0.3*std(acc)to each original acceleration signal
NOISYsignal_sub=NOISYsignal-mean(NOISYsignal);
ORIGINALacc_sub=acc(indx,:)-m_acc;
[PSD,Freq]=pwelch(ORIGINALacc_sub,1000,[],F,Fs); %PSD of original signal
[PSD1,Freq1]= pwelch(NOISYsignal_sub,1000,[],F,Fs);%PSD of noisy+original signal
PSD_acc=[PSD_acc,PSD];
PSD_acc1=[PSD_acc1,PSD1];
Freq_acc=[Freq_acc,Freq];
Freq_acc1=[Freq_acc1,Freq1];
Cutoff = 20;
Indx_signal = find(Freq1<=Cutoff);
Indx_noise = find(Freq1>Cutoff);
PSD_signal1 = sum(PSD1(Indx_signal));
PSD_noise1 = sum(PSD1(Indx_noise));
SNR1 = PSD_signal1/PSD_noise1;
%True SNR (this is the code we will use to find true SNR
% true_noise = signal-sum_harmonics;
% var_signal_true = var(sum_harmonics);
% var_noise_true = var(signal-sum_harmonics);
% SNR_true = var_signal_true/var_noise_true;
%moving average
for m=1:length(NOISYsignal_sub)-Win
smooth_acc(m)=mean(NOISYsignal_sub(m:m+Win));
win_time(m)=time(indx(m+win_center));
end
end
if Fig==1
plot(time(indx),NOISYsignal_sub,'b-');hold on %plotting of noisy signal: by taking minus mean--> zero
plot(time(indx),acc(indx),'r-',"LineWidth",1); %plotting of original data: there is an off set in the original data, it happens due to incorrect offset
plot(win_time,smooth_acc,'g',"LineWidth",2); %plotting of moving average in green
xlabel('time');ylabel('acceleration');
title('Noisy acceleration with moving average')
for n=1:N
subplot(2,2,n)
plot(Freq_acc(:,n),PSD_acc(:,n));hold on
plot(Freq_acc1(:,n),PSD_acc1(:,n));
title(Labels_acc{n});
xlabel('Frequency(Hz)');ylabel('PSD');
legend('Original Signal','Noisy signal');
end
end
end
hanna Eid on 22 May 2022
How do you make PSD frequency into a signal?

Abhimenyu on 8 Dec 2023
Hi Hanna,
I understand that you want to find the true Signal-to-noise ratio of the acceleration signal using the sum of harmonics. The “harmonic” function of MATLAB can be used to find the sum of harmonics. The “harmonic” function is defined as the sum of reciprocals of the first “n” positive integers where “n” is the input.
To find the sum of harmonics of a signal, first MATLAB’sfft function can be used to compute the discrete Fourier transform of the signal, and then the “harmonic function can be used to sum the reciprocals of the magnitudes of the frequency components. Please refer to the example code below to understand more:
% Compute the discrete Fourier transform of the signal
X = fft(signal);
% Get the magnitudes of the frequency components
Mag = abs(X);
% Sum the reciprocals of the magnitudes using the "harmonic" function
sum_harmonics = harmonic(M);
The variable “sum_harmonics will contain the sum of harmonics of the signal. This value can be used to compute the true signal-to-noise ratio (SNR) of the signal according to the code mentioned by you.
Please refer to the below given MATLAB documentation links to explore more on the “harmonic” and the “fft function respectively:
I hope this helps to resolve the query.
Thanks,
Abhimenyu