www.gusucode.com > signal 案例源码程序 matlab代码 > signal/FilterIntroductionExample.m
%% Practical Introduction to Digital Filtering % % This example shows how to design, analyze, and apply a digital filter to % your data. It will help you answer questions such as: how do I compensate % for the delay introduced by a filter?, How do I avoid distorting my % signal?, How do I remove unwanted content from my signal?, How do I % differentiate my signal?, and How do I integrate my signal? % % Filters can be used to shape the signal spectrum in a desired way or to % perform mathematical operations such as differentiation and integration. % In what follows you will learn some practical concepts that will ease the % use of filters when you need them. % % This example focuses on applications of digital filters rather than on % their design. If you want to learn more about how to design digital % filters see the % <../../../signal/sigdemos/html/FilterDesignIntroductionExample.html % "Practical Introduction to Digital Filter Design"> example. % Copyright 2012-2013 The MathWorks, Inc. %% Compensating for Delay Introduced by Filtering % % Digital filters introduce delay in your signal. Depending on the filter % characteristics, the delay can be constant over all frequencies, or it % can vary with frequency. The type of delay determines the actions you % have to take to compensate for it. The *grpdelay* function allows you to % look at the filter delay as a function of frequency. Looking at the % output of this function allows you to identify if the delay of the filter % is constant or if it varies with frequency (i.e. if it is % frequency-dependent). % % Filter delay that is constant over all frequencies can be easily % compensated for by shifting the signal in time. FIR filters usually have % constant delay. On the other hand, delay that varies with frequency % causes phase distortion and can alter a signal waveform significantly. % Compensating for frequency-dependent delay is not as trivial as for the % constant delay case. IIR filters introduce frequency-dependent delay. % % *Compensating for Constant Filter Delay* % % As mentioned before, you can measure the group of delay of the filter to % verify that it is a constant function of frequency. You can use the % *grpdelay* function to measure the filter delay, D, and compensate for % this delay by appending D zeros to the input signal and shifting the % output signal in time by D samples. % % Consider a noisy electrocardiogram signal that you want to filter to % remove high frequency noise above 75 Hz. You want to apply an FIR lowpass % filter and compensate for the filter delay so that the noisy and filtered % signals are aligned correctly and can be plotted on top of each other for % comparison. Fs = 500; % sample rate in Hz N = 500; % number of signal samples rng default; x = ecg(N)'+0.25*randn(N,1); % noisy waveform t = (0:N-1)/Fs; % time vector % Design a 70th order lowpass FIR filter with cutoff frequency of 75 Hz. Fnorm = 75/(Fs/2); % Normalized frequency df = designfilt('lowpassfir','FilterOrder',70,'CutoffFrequency',Fnorm); %% % Plot the group delay of the filter to verify that it is constant across % all frequencies indicating that the filter is linear phase. Use the group % delay to measure the delay of the filter. grpdelay(df,2048,Fs) % plot group delay D = mean(grpdelay(df)) % filter delay in samples %% % Before filtering, append D zeros at the end of the input data vector, x. % This ensures that all the useful samples are flushed out of the filter, % and that the input signal and the delay-compensated output signal have % the same length. Filter the data and compensate for the delay by shifting % the output signal by D samples. This last step effectively removes the % filter transient. y = filter(df,[x; zeros(D,1)]); % Append D zeros to the input data y = y(D+1:end); % Shift data to compensate for delay figure plot(t,x,t,y,'r','linewidth',1.5); title('Filtered Waveforms'); xlabel('Time (s)') legend('Original Noisy Signal','Filtered Signal'); grid on axis tight %% % *Compensating for Frequency-Dependent Delay* % % Frequency-dependent delay causes phase distortion in the signal. % Compensating for this type of delay is not as trivial as for the constant % delay case. If your application allows off-line processing, you can % remove the frequency-dependent delay by implementing zero-phase filtering % using the *filtfilt* function. *filtfilt* performs zero-phase filtering % by processing the input data in both the forward and reverse directions. % The main effect is that you obtain zero-phase distortion, i.e., you % filter data with an equivalent filter that has a constant delay of 0 % samples. Other effects are that you get a filter transfer function which % equals the squared magnitude of the original filter transfer function, % and a filter order that is double the order of the original filter. % % Consider the ECG signal defined in the previous section. Filter this % signal with and without delay compensation. % Design a 7th order lowpass IIR elliptic filter with cutoff frequency % of 75 Hz. Fnorm = 75/(Fs/2); % Normalized frequency df = designfilt('lowpassiir',... 'PassbandFrequency',Fnorm,... 'FilterOrder',7,... 'PassbandRipple',1,... 'StopbandAttenuation',60); %% % Plot the group delay of the filter and notice that it varies with % frequency indicating that the filter delay is frequency-dependent. grpdelay(df,2048,'half',Fs) %% % Filter the data and look at the effects of each filter implementation on % the time signal. y1 = filter(df,x); % non-linear phase filter - no delay compensation y2 = filtfilt(df,x); % zero-phase implementation - delay compensation figure plot(t,x); hold on plot(t,y1,'r','linewidth',1.5); plot(t,y2,'g','linewidth',1.5); title('Filtered Waveforms'); xlabel('Time (s)') legend('Original Signal','Non-linear phase IIR output',... 'Zero-phase IIR output'); ax = axis; axis([0.25 0.55 ax(3:4)]) grid on %% % Notice how zero-phase filtering effectively removes the filter delay. % % Zero-phase filtering is a great tool if your application allows for the % non-causal forward/backward filtering operations, and for the change of % the filter response to the square of the original response. % % Filters that introduce constant delay are linear phase filters. Filters % that introduce frequency-dependent delay are non-linear phase filters. %% Removing Unwanted Spectral Content from a Signal % % Filters are commonly used to remove unwanted spectral content from a % signal. You can choose from a variety of filters to do this. You choose a % lowpass filter when you want to remove high frequency content, or a % highpass filter when you want to remove low frequency content. You can % also choose a bandpass filter to remove low and high frequency content % while leaving an intermediate band of frequencies intact. You choose a % bandstop filter when you want to remove frequencies over a given band. % % Consider an audio signal that has a power-line hum and white noise. The % power-line hum is caused by a 60 Hz tone. White noise is a signal that % exists across all the audio bandwidth. % % Load the audio signal. Fs = 44100; % Sample rate y = audioread('noisymusic.wav'); %% % Plot the power spectrum of the signal. The red triangular marker shows % the strong 60 Hz tone interfering with the audio signal. [P,F] = pwelch(y,ones(8192,1),8192/2,8192,Fs,'power'); helperFilterIntroductionPlot1(F,P,[60 60],[-9.365 -9.365],... {'Original signal power spectrum', '60 Hz Tone'}) %% % You can first remove as much white noise spectral content as possible % using a lowpass filter. The passband of the filter should be set to a % value that offers a good trade-off between noise reduction and audio % degradation due to loss of high frequency content. Applying the lowpass % filter before removing the 60 Hz hum is very convenient since you will be % able to downsample the band-limited signal. The lower rate signal will % allow you to design a sharper and narrower 60 Hz bandstop filter with a % smaller filter order. % % Design a lowpass filter with passband frequency of 1 kHz, and stopband % frequency of 1.4 kHz. Choose a minimum order design. Fp = 1e3; % Passband frequency in Hz Fst = 1.4e3; % Stopband frequency in Hz Ap = 1; % Passband ripple in dB Ast = 95; % Stopband attenuation in dB % Design the filter df = designfilt('lowpassfir','PassbandFrequency',Fp,... 'StopbandFrequency',Fst,'PassbandRipple',Ap,... 'StopbandAttenuation',Ast,'SampleRate',Fs); % Analyze the filter response hfvt = fvtool(df,'Fs',Fs,'FrequencyScale','log',... 'FrequencyRange','Specify freq. vector','FrequencyVector',F); %% % Filter the data and compensate for delay D = mean(grpdelay(df)); % filter delay ylp = filter(df,[y; zeros(D,1)]); ylp = ylp(D+1:end); close(hfvt) %% % Look at the spectrum of the lowpass filtered signal. Note how the % frequency content above 1400 Hz has been removed. [Plp,Flp] = pwelch(ylp,ones(8192,1),8192/2,8192,Fs,'power'); helperFilterIntroductionPlot1(F,P,Flp,Plp,... {'Original signal','Lowpass filtered signal'}) %% % From the power spectrum plot above, you can see that the maximum % non-negligible frequency content of the lowpass filtered signal is at % 1400 Hz. By the sampling theorem, a sample frequency of 2*1400 = 2800 Hz % would suffice to represent the signal correctly, you however, are using a % sample rate of 44100 Hz which is a waste since you will need to process % more samples than those necessary. You can downsample the signal to % reduce the sample rate and reduce the computational load by reducing the % number of samples that you need to process. A lower sample rate will also % allow you to design a sharper and narrower bandstop filter, needed to % remove the 60 Hz noise, with a smaller filter order. % % Downsample the lowpass filtered signal by a factor of 10 to obtain a % sample rate of Fs/10 = 4.41 kHz. Plot the spectrum of the signal before % and after downsampling. Fs = Fs/10; yds = downsample(ylp,10); [Pds,Fds] = pwelch(yds,ones(8192,1),8192/2,8192,Fs,'power'); helperFilterIntroductionPlot1(F,P,Fds,Pds,... {'Signal sampled at 44100 Hz', 'Downsampled signal, Fs = 4410 Hz'}) %% % Now remove the 60 Hz tone using an IIR bandstop filter. Let the stopband % have a width of 4 Hz centered at 60 Hz. We choose an IIR filter to % achieve a sharp frequency notch, small passband ripple, and a relatively % low order. Process the data using *filtfilt* to avoid phase distortion. % Design the filter df = designfilt('bandstopiir','PassbandFrequency1',55,... 'StopbandFrequency1',58,'StopbandFrequency2',62,... 'PassbandFrequency2',65,'PassbandRipple1',1,... 'StopbandAttenuation',60,'PassbandRipple2',1,... 'SampleRate',Fs,'DesignMethod','ellip'); % Analyze the magnitude response hfvt = fvtool(df,'Fs',Fs,'FrequencyScale','log',... 'FrequencyRange','Specify freq. vector','FrequencyVector',Fds(Fds>F(2))); %% % Perform zero-phase filtering to avoid distortion. ybs = filtfilt(df,yds); %% % Finally, upsample the signal to bring it back to the original audio % sample rate of 44.1 kHz which is compatible with audio soundcards. yf = interp(ybs,10); Fs = Fs*10; %% % Take a final look at the spectrum of the original and processed signals. % Notice how the high frequency noise floor and the 60 Hz tone have been % attenuated by the filters. [Pfinal,Ffinal] = pwelch(yf,ones(8192,1),8192/2,8192,Fs,'power'); close(hfvt) helperFilterIntroductionPlot1(F,P,Ffinal,Pfinal,... {'Original signal','Final filtered signal'}) %% % Listen to the signal before and after processing. As mentioned above, the % end result is that you have effectively attenuated the 60 Hz hum and the % high frequency noise on the audio file. % Play the original signal hplayer = audioplayer(y, Fs); play(hplayer); % Play the noise-reduced signal hplayer = audioplayer(yf, Fs); play(hplayer); %% Differentiating a Signal % % The MATLAB *diff* function differentiates a signal with the drawback that % you can potentially increase the noise levels at the output. A better % option is to use a differentiator filter that acts as a differentiator in % the band of interest, and as an attenuator at all other frequencies, % effectively removing high frequency noise. % % As an example, analyze the speed of displacement of a building floor % during an earthquake. Displacement or drift measurements were recorded on % the first floor of a three story test structure under earthquake % conditions and saved in the quakedrift.mat file. The length % of the data vector is 10e3, the sample rate is 1 kHz, and the units of % the measurements are cm. % % Differentiate the displacement data to obtain estimates of the speed and % acceleration of the building floor during the earthquake. Compare the % results using diff and an FIR differentiator filter. load quakedrift.mat Fs = 1000; % sample rate dt = 1/Fs; % time differential t = (0:length(drift)-1)*dt; % time vector %% % Design a 50th order differentiator filter with a passband frequency of % 100 Hz which is the bandwidth over which most of the signal energy is % found. Set the stopband frequency of the filter to 120 Hz. df = designfilt('differentiatorfir','FilterOrder',50,... 'PassbandFrequency',100,'StopbandFrequency',120,... 'SampleRate',Fs); %% % The *diff* function can be seen as a first order FIR filter with response % $H(Z) = 1 - Z^{-1}$. Use FVTool to compare the magnitude response of the % 50th order differentiator FIR filter and the response of the *diff* % function. Clearly, both responses are equivalent in the passband region % (from 0 to 100 Hz). However, in the stopband region, the 50th order % filter attenuates components while the diff response amplifies % components. This effectively increases the levels of high frequency % noise. hfvt = fvtool(df,[1 -1],1,'MagnitudeDisplay','zero-phase','Fs',Fs); legend(hfvt,'50th order FIR differentiator','Response of diff function'); %% % Differentiate using the *diff* function. Add zeros to compensate for the % missing samples due to the diff operation. v1 = diff(drift)/dt; a1 = diff(v1)/dt; v1 = [0; v1]; a1 = [0; 0; a1]; %% % Differentiate using the 50th order FIR filter and compensate for delay. D = mean(grpdelay(df)); % filter delay v2 = filter(df,[drift; zeros(D,1)]); v2 = v2(D+1:end); a2 = filter(df,[v2; zeros(D,1)]); a2 = a2(D+1:end); v2 = v2/dt; a2 = a2/dt^2; %% % Plot a few data points of the floor displacement. Plot also a few data % points of the speed and acceleration as computed with diff and with the % 50th order FIR filter. Notice how the noise has been slightly amplified % in the speed estimates and largely amplified in the acceleration % estimates obtained with *diff*. helperFilterIntroductionPlot2(t,drift,v1,v2,a1,a2) %% Integrating a Signal % % A leaky integrator filter is an all-pole filter with transfer function % $H(Z) = 1/[1-cZ^{-1}]$ where $c$ is a constant that must be smaller than % 1 to ensure stability of the filter. It is no surprise that as $c$ % approaches one, the leaky integrator approaches the inverse of the *diff* % transfer function. Apply the leaky integrator to the acceleration and % speed estimates obtained in the previous section to get back the speed % and the drift respectively. Use the estimates obtained with the *diff* % function since they are noisier. % % Use a leaky integrator with $a = 0.999$. Plot the magnitude response of % the leaky integrator filter. Notice that the filter acts as a lowpass % filter effectively eliminating high frequency noise. close(hfvt) fvtool(1,[1 -.999],'Fs',Fs) %% % Filter the velocity and acceleration with the leaky integrator. v_original = v1; a_original = a1; d_leakyint = filter(1,[1 -0.999],v_original); v_leakyint = filter(1,[1 -0.999],a_original); % Multiply by time differential d_leakyint = d_leakyint * dt; v_leakyint = v_leakyint * dt; %% % Plot the displacement and speed estimates and compare to the original % signals v1 and a1. helperFilterIntroductionPlot3(t,drift,d_leakyint,v_original,v_leakyint) %% % You can also integrate a signal using the *cumsum* and *cumtrapz* % functions. Results will be similar to those obtained with the leaky % integrator. %% Conclusions % % In this example you learned about linear and nonlinear phase filters and % you learned how to compensate for the phase delay introduced by each % filter type. You also learned how to apply filters to remove unwanted % frequency components from a signal, and how to downsample a signal after % limiting its bandwidth with a lowpass filter. Finally, you learned how to % differentiate and integrate a signal using digital filter designs. % Throughout the example you also learned how to use analysis tools to look % at the response and group delay of your filters. %% Further Reading % % For more information on filter applications see the Signal Processing % Toolbox. For more information on how to design digital filters see the % <../../../signal/sigdemos/html/FilterDesignIntroductionExample.html % "Practical Introduction to Digital Filter Design"> example. % % References: % J.G. Proakis and D. G. Manolakis, "Digital Signal Processing. % Principles, Algorithms, and Applications", Prentice-Hall, 1996. % % S. J. Orfanidis, "Introduction To Signal Processing", Prentice-Hall, % 1996. %% Appendix % The following helper functions are used in this example. % % * <matlab:edit('helperFilterIntroductionPlot1.m') helperFilterIntroductionPlot1.m> % * <matlab:edit('helperFilterIntroductionPlot2.m') helperFilterIntroductionPlot2.m> % * <matlab:edit('helperFilterIntroductionPlot3.m') helperFilterIntroductionPlot3.m>