The noisy signal is processed by mean filtering

Using matlab to filter the noisy sine wave signal
This is my little job. There may be all kinds of small errors in the program. Firstly, the sinusoidal signal is generated, different signal-to-noise ratios are set, and the moving average filter is designed to filter the noisy signal.
For a fixed signal-to-noise ratio, the influence of the order of the moving average filter on the filtering performance is analyzed and discussed, the root mean square error of the denoised signal and the real signal is calculated, the relationship curve between the root mean square error and the order of the filter is drawn, and the comparison curve between the filtered signal and the original signal and the noisy signal under the optimal filter order is drawn.
Principle:
Firstly, generate a sine wave signal, set corresponding parameters, set each parameter as a variable, superimpose noise on the signal, superimpose noise on the signal using awgn function, set the range of signal-to-noise ratio, and define corresponding variables to facilitate parameter change in the later stage.
Calculate and store the matrix containing noise signal through the above parameters, select a fixed signal-to-noise ratio, filter this group of data, set the order of different mean filter and the order storage matrix of filter, and set the storage matrix of filtered data to filter the data
Finally, the filter function is designed to filter the signal with filters of different orders, and the error calculation function is designed to calculate the minimum mean square error under different filter orders
code:

clear,clc,close all;
f=50;%frequency
fs=10000;%sampling frequency 
Ts=1/fs;%Sampling interval
N=5000;%Sampling points
n=1:N;
y=sin(2*pi*f*n*Ts);

SNR_strat=-10;
SNR_end=50;
SNR_step=1;
SNR=SNR_strat:SNR_step:SNR_end;
len_add_noise_row=N;
len_add_noise_col=length(SNR)+1;
data_noise=zeros(len_add_noise_row,len_add_noise_col);
data_noise(:,1)=y;
figure,plot(data_noise(:,1));title('sin(2*pi*f*n*Ts)');
cnt=2;
for loop=SNR
    data_noise(:,cnt)=awgn(data_noise(:,1),loop,'measured');
    figure,plot(data_noise(:,cnt));title(strcat('SNR=',num2str(loop)));
    cnt=cnt+1;
end
close all;

SNR_sle_idex=37;
%SNR_sle_idex=31;
SNR_sle=SNR(SNR_sle_idex);
data_sle=data_noise(:,SNR_sle_idex+1);


M_start=2;%The minimum starting value is 2
M_step=1;
M_end=80;
M=M_start:M_step:M_end;


data_filter=zeros(N,length(M)+2);
data_filter(:,1)=y;

figure,plot(data_noise(:,1));title('sin(2*pi*f*n*Ts)');
figure,plot(data_sle);title(strcat('Filtered signal waveform',strcat('The signal-to-noise ratio is',num2str(SNR_sle))));
data_filter(:,2:end)=filter_M_mean(SNR_sle,data_sle ,M);

error_data_filter_matrix=zeros(1,length(M)+1);
error_data_filter_matrix=error_M_data(data_filter);
M_add_y=[ 1 M];

min_data=min(error_data_filter_matrix);
index_min=find(error_data_filter_matrix==min_data);
close all;
figure,plot(M_add_y,error_data_filter_matrix);
xlabel('M');ylabel('Mean square error');
title(strcat(strcat(strcat(strcat(strcat('When the signal-to-noise ratio is',num2str(SNR_sle)),'dB When the order of the filter is M='),num2str(M(index_min-1))),'The root mean square error is the smallest, and the root mean square error is'),num2str(min_data)));
text(M(index_min-1),min_data,strcat('M=',num2str(M(index_min-1))))
hold on;
plot(M(index_min-1),min_data,'r+');
disp(strcat(strcat(strcat(strcat(strcat('When the signal-to-noise ratio is',num2str(SNR_sle)),'dB When the order of the filter is M='),num2str(M(index_min-1))),'The root mean square error is the smallest, and the root mean square error is'),num2str(min_data)));

figure,subplot(311);plot(data_filter(:,1));xlabel('t');ylabel('range');title('original signal ');grid on;
subplot(312);plot(data_filter(:,2));xlabel('t');ylabel('range');title(strcat(strcat('The signal-to-noise ratio is',num2str(SNR_sle)),'dB In the case of noisy signal diagram'));
axis([-inf inf -inf inf]);grid on;
subplot(313);plot(data_filter(:,index_min+1))
xlabel('t');ylabel('range');title(strcat(strcat(strcat('Filter order M=',num2str(M(index_min-1))),'The minimum error of mean square value is'),num2str(min_data)));
axis([-inf inf -inf inf]);grid on;

M_optimal=M(index_min-1);
SNR_2=SNR;
data_noise_2=data_noise;
data_noise_3=M_optimal_filter_data(data_noise_2,M_optimal);
%data_noise_3  The first column is the original data, followed by noisy signals%Matrix for storing mean square error
data_M_optimal_all_SNR=zeros(1,length(SNR_2));
data_M_optimal_all_SNR=error_M_SNR(data_noise_3);
figure,plot(SNR_2,data_M_optimal_all_SNR);xlabel('SNR');ylabel('Mean square error');title(strcat(strcat('Mean filter order M=',num2str(M_optimal)),'Variation curve of bit error rate with signal-to-noise ratio'));

The rest of the called functions

function [ data_noise_3 ] = M_optimal_filter_data( data_noise_2,M_optimal )
[M,N]=size(data_noise_2);
data_noise_3=zeros(M,N)+inf;
data_noise_3(:,1)=data_noise_2(:,1);
for j=2:N
    for i=1:M
        if(i==1|i>(N-(M_optimal-2)))
            data_noise_3(i,j)=data_noise_2(i,j);%Take the original value directly for the boundary value
        end
        if(~(i==1|i>(N-(M_optimal-2))))
            data_noise_3(i,j)=0;
            for k=-1:M_optimal-2
            data_noise_3(i,j)=data_noise_3(i,j)+data_noise_2(i+k,j);
            end     
            data_noise_3(i,j)=data_noise_3(i,j)/M_optimal;
        end
    end
end
end

function [ data_return ] = filter_M_mean(SNR_sle,data_sle ,M )
data_return=zeros(length(data_sle),length(M)+1)+inf;
data_return(:,1)=data_sle;
M_start=M(1);
M_end=M(end);
M_step=M(end)-M(end-1);
cnt=2;
for loop=M_start:M_step:M_end
    for i=1:length(data_sle)
    if(i==1|i>(length(data_sle)-(loop-2)))
        data_return(i,cnt)=(data_return(i,1))/loop;
    else
        data_return(i,cnt)=0;
        for j=-1:loop-2
            data_return(i,cnt)=data_return(i,cnt)+data_return(i+j,1);
        end
        data_return(i,cnt)=data_return(i,cnt)/loop;
    end
    end
    figure,plot(data_return(:,cnt));title(strcat(strcat('SNR=',num2str(SNR_sle)),strcat('Filter order',num2str(M(cnt-1)))));
    cnt=cnt+1;
end
disp(data_return);
end
function [ data_M_optimal_all_SNR ] = error_M_SNR( data_noise_3 )
[M,N]=size(data_noise_3);
data_M_optimal_all_SNR=zeros(1,N-1);
for j=2:N
    sum=0;
    for i=1:M
        sum=sum+(data_noise_3(i,j)-data_noise_3(i,1))^2;
    end
    data_M_optimal_all_SNR(1,j-1)=sqrt(sum/M);
end
end
function [ error_data_filter_matrix ] = error_M_data( data_filter )
[M,N]=size(data_filter);
error_data_filter_matrix=zeros(1,N-1)+inf;
for i=2:N
    sum=0;
    for j=1:M
        sum=sum+(data_filter(j,i)-data_filter(j,1))^2;
    end
    error_data_filter_matrix(1,i-1)=sqrt(sum/M);
end
end


Given the order of filter, change the signal-to-noise ratio of noisy signal, analyze the filtering performance under different signal-to-noise ratio, and draw the variation curve of root mean square error with signal-to-noise ratio. (see above for code)

The second problem is: given a flow signal, find the period of the signal
code:

%Read data
close all;
[num,txt,raw]=xlsread('01.xlsx');
[M,N]=size(num);
figure,plot(num(:,1));xlabel('t');ylabel('flow');title('Signal diagram');
figure,plot(num(:,2));xlabel('t');ylabel('speed');title('Velocity diagram');
disp('Filter the flow signal');
M_start=2;
M_step=5;
M_end=200;
M=M_start:M_step:M_end;
[data_filter_traffic]=main(num(:,1)',M);
index_T=19;%This variable must be greater than or equal to 2
M_T=M(index_T-1);
disp(strcat(strcat('Select the filter order as',num2str(M_T)),'Calculation period of operation results of first-order mean filter'));
close all;
figure,plot(num(:,1));xlabel('t');ylabel('flow');title('Signal diagram');
figure,plot(data_filter_traffic(:,index_T));xlabel('t');ylabel('flow');title(strcat(strcat('Select the filter order as',num2str(M_T)),'Filtering result diagram'));grid on; 
data_T=data_filter_traffic(:,index_T);
[T1,T_traffic,data_T1,T_traffic1]=Find_T(data_T);
disp('The maximum value detection method has high requirements for the filter. This code requires the average filter order to be above 87');
disp(strcat('Calculation cycle of maximum detection value:',num2str(T1)));
T2=Find_T_other(data_T,M_T);
disp(strcat('Zero crossing detection calculation cycle:',num2str(T2)));
function [data_filter] = main( y ,M)
N=length(y);
data_sle=y;
data_filter=zeros(N,length(M)+1);
data_filter(:,1)=y';
data_filter(:,1:end)=filter_M_mean(y ,M);
end
function [ T_traffic1 ,T_return] = sub_distance( T_traffic,data_T)
data_T=data_T';
N=length(T_traffic);
T_traffic1=zeros(1,N-1)+inf;
for loop=2:N
    T_traffic1(1,loop-1)=T_traffic(loop)-T_traffic(loop-1);
end
data_find=find(T_traffic1>mean(T_traffic1));
data1=data_T';
figure,plot(data_T);grid on;xlabel('t');ylabel('flow');hold on;title('Results after classification');hold on;
plot(T_traffic(data_find),data_T(T_traffic(data_find)),'r*');
hold on;
plot(T_traffic(end),data_T(T_traffic(end)),'r*');
T_return=round(length(data_T)/(length(data_find)+1));
end
function [ data_noise_3 ] = M_optimal_filter_data( data_noise_2,M_optimal )
[M,N]=size(data_noise_2);
data_noise_3=zeros(M,N)+inf;
data_noise_3(:,1)=data_noise_2(:,1);
for j=2:N
    for i=1:M
        if(i==1|i>(N-(M_optimal-2)))
            data_noise_3(i,j)=data_noise_2(i,j);
        end
        if(~(i==1|i>(N-(M_optimal-2))))
            data_noise_3(i,j)=0;
            for k=-1:M_optimal-2
            data_noise_3(i,j)=data_noise_3(i,j)+data_noise_2(i+k,j);
            end     
            data_noise_3(i,j)=data_noise_3(i,j)/M_optimal;
        end
    end
end
end
function [ T ] = Find_T_other( data_T,M_T )
N=length(data_T);
data_T_sub_mean=data_T-mean(data_T);
data_T_sub_mean=data_T_sub_mean';
cnt=0;
index=inf;
for loop=2:N-1
    if(((data_T_sub_mean(1,loop-1)<0)&&(data_T_sub_mean(1,loop+1)>0))||((data_T_sub_mean(1,loop-1)>0)&&(data_T_sub_mean(1,loop+1)<0)))
        cnt=cnt+1;
        index=[index loop];
    end
end
index_1=index(:,2:end);
data=data_T';
data=data(index_1);
figure,plot(data_T_sub_mean);xlabel('t');ylabel('flow');
title(strcat(strcat(strcat('Select the filter order as',num2str(M_T)),'Filtering result diagram'),'The result of subtracting the mean value of the filtered signal,gules+Represents the zero point determined by the system'));grid on;
hold on;
plot(index_1,data-mean(data_T),'r+');
% disp(cnt);
% disp(index_1);
 T=round(N/(cnt/4));
end

function [ T_return,T_traffic,data_T1,T_traffic1] = Find_T( data_T )
M=length(data_T);
data_T=data_T';
T_traffic=-1;
M_delay=10;
for loop=2:M-1
    flag=1;
    if(((data_T(loop)>data_T(loop-1))&&(data_T(loop)>data_T(loop+1))&&(data_T(loop)>mean(data_T))))
        T_traffic=[T_traffic loop];
        flag=0;
    end
    if (data_T(loop)>mean(data_T)&&(flag&&((loop+M_delay<M-1)&&(loop-M_delay>2)&&(((data_T(loop)>data_T(loop-1))&&(data_T(loop)>data_T(loop+1))&&(data_T(loop)>mean(data_T)))||(data_T(loop)>data_T(loop-M_delay))&&(data_T(loop)>data_T(loop+M_delay))))))
       T_traffic=[T_traffic loop];
    end
end
T_traffic=T_traffic(2:end);
N=length(T_traffic);
data=data_T';
data=data(T_traffic);
T_t=zeros(1,N-1)-1;
for i=1:N-1
    T_t(1,i)=T_traffic(1,i+1)-T_traffic(1,i);
end
figure,plot(data_T);grid on;xlabel('t');ylabel('flow');hold on;title('gules+Represents the maximum value of an interval detected by the system');
plot(T_traffic,data,'r+');hold on;
data_T1=data_T(T_traffic);
[T_traffic1,T_return]=sub_distance(T_traffic,data_T);%T_traffic1 Is the difference matrix
end


function [ data_return ] = filter_M_mean(data_sle ,M )
data_return=zeros(length(data_sle),length(M)+1)+inf;
data_return(:,1)=data_sle';
M_start=M(1);
M_end=M(end);
M_step=M(end)-M(end-1);
cnt=2;
for loop=M_start:M_step:M_end
    for i=1:length(data_sle)
    if(i==1|i>(length(data_sle)-(loop-2)))
        data_return(i,cnt)=(data_return(i,1))/loop;
    else
        data_return(i,cnt)=0;
        for j=-1:loop-2
            data_return(i,cnt)=data_return(i,cnt)+data_return(i+j,1);
        end
        data_return(i,cnt)=data_return(i,cnt)/loop;
    end
    end
    %figure,plot(data_return(:,cnt));title(strcat('Filter order',num2str(M(cnt-1))));
    cnt=cnt+1;
end
end


function [ data_M_optimal_all_SNR ] = error_M_SNR( data_noise_3 )
[M,N]=size(data_noise_3);
data_M_optimal_all_SNR=zeros(1,N-1);
for j=2:N
    sum=0;
    for i=1:M
        sum=sum+(data_noise_3(i,j)-data_noise_3(i,1))^2;
    end
    data_M_optimal_all_SNR(1,j-1)=sqrt(sum/M);
end
end


function [ error_data_filter_matrix ] = error_M_data( data_filter )
[M,N]=size(data_filter);
error_data_filter_matrix=zeros(1,N-1)+inf;
for i=2:N
    sum=0;
    for j=1:M
        sum=sum+(data_filter(j,i)-data_filter(j,1))^2;
    end
    error_data_filter_matrix(1,i-1)=sqrt(sum/M);
end
end

Keywords: MATLAB

Added by BZorch on Sun, 24 Oct 2021 14:36:42 +0300