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

2 次查看（过去 30 天）
hanna Eid
hanna Eid 2022-5-22
回答： Abhimenyu 2023-12-8
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");
load Index_table.mat
Fig=1
Fs=100; %sampling frequency
F= Fs/2;
for i=1:length(files)
load(files(i).name);%task 2
indx=find(time>=Index_tabel2{i,1}& time<=Index_tabel2{i,2}); %task 3
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
figure;%task 8: moving average
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')
figure;%task 5b
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

请先登录，再回答此问题。

回答（1 个）

Abhimenyu
Abhimenyu 2023-12-8
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’s “fft” 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

类别

在 Help Center 和 File Exchange 中查找有关 Discrete Fourier and Cosine Transforms 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by