Lynn's integer filter has a linearly increasing trend I have a PPG signal sampled at 50Hz and I want to filter the signal using Lynn's integer filter. I.e

\begin{equation}
  H(z)=\frac{(1-1/z^m)^2}{(1-1/z)^2}
\end{equation}

Cutoff frequency is 5Hz, so this gives value of $m \approx 3 $ (it's a 3db cutoff) and hence the difference equation is

\begin{equation}
y[n]-2y[n-1]+y[n-2]=x[n]-2x[n-3]+x[n-6]
\end{equation}

Now applying this filter to the signal, my signal is filtered but there is some linear component that's making baseline to increase linearly.

Any suggestion where I am going wrong?

Data is processed on cortex M0 based bluetooth module and transmitted to phone. I am plotting csv file in Matlab.

Code of cortex M0:

float process_ppg_iir(uint16_t raw_data_ppg)
{
    int b0=1, b1=-2, b2=1, a1=-2;
    int a2 = 1;

    if( global_sample_counter <= 7-1 )
    {
        ppg_raw[global_sample_counter] = raw_data_ppg;
        if(global_sample_counter <=3-1)
        ppg_processed[global_sample_counter] = raw_data_ppg;
        global_sample_counter++;
    }
    else
    {
        // sfit values to left

        int i=0;

        for(i=0 ;i <=5;i++)
        {
            ppg_raw[i] = ppg_raw[i+1];
        }

        ppg_processed[0] = ppg_processed[1];
        ppg_processed[1] = ppg_processed[2];

        ppg_raw[6] = raw_data_ppg;
        ppg_processed[2] = (b0*ppg_raw[6]+b1*ppg_raw[3]+b2*ppg_raw[0])- a1*ppg_processed[1]-a2*ppg_processed[0];

    }
        return ppg_processed[2];
}

Your lowpass filter is the cascade of two 3-point moving average filters each implemented as recursive running sum filters. (It's gain at 5 Hz is -2.4 dB.) I can't interpret your code but I'm guessing your feedback operations are incorrect. To test your code, make sure all accumulator variables are set to zero. Then apply a unit impulse to your filter and see if its impulse response is [1, 2, 3, 2, 1, 0, 0, 0, ...] as it should be.  If you change your $H(z)$ exponents from $2$ to $1$, then your filter's impulse response should be [1, 1, 1, 0, 0,0,  ...]. If your impulse responses do not agree with the above, then finding out why they don't will solve your problem.


How to plot magnitude and phase response of 2 cascaded filters in Matlab? How would I go about plotting a magnitude and phase response of a system that consists of two cascaded 2nd order Butterworth filters in Matlab? Filters are the same.

% Filter coeffs:
[B, A] = butter( 2, 1000/8000*2, 'high' ); % Cutoff at 1000 Hz, fs 8000 Hz.

% Cascade 2 filters: input -> filter -> filter -> output :
[out1, state1 ] = filter( B, A, in,   state1);
[out2, state2 ] = filter( B, A, out1, state2);

Using MATLAB/Octave as the tool, the following approach lets you plot the magnitude & phase samples of the DTFT of the cascade of the two discrete-time LTI filters using their LCCDE coefficient vectors $b[k]$, and $a[k]$, assuming they have LCCDE representations..

Since the cascade LTI system is described as:
$$h[n] = h_1[n] \star h_2[n]$$ and consequently
$$ H(z) = H_1(z) H_2(z) $$

Where the $h[n]$ is the impulse response of the cascade system and $H(z)$ is its transfer function (Z-transform of impulse response). Those $h_1[n]$ and $h_2[n]$ refer to the impulse responses of the individual systems that make up the cascade.

Assuming that individual systems do have rational Z-transforms, then we have the following expression for the composite (cascade) system Z-transfom:
$$
\begin{align}
H(z)  &=\frac{ \sum_{k=0}^{M} B[k]z^{-k} } {\sum_{k=0}^{N} A[k]z^{-k}} = H_1(z) H_2(z) \\ \\
&=
\left(  \frac{ \sum_{k=0}^{M_1} b_1[k]z^{-k} } {\sum_{k=0}^{N_1} a_1[k]z^{-k}} \right)
\left( \frac{ \sum_{k=0}^{M_2} b_2[k]z^{-k} } {\sum_{k=0}^{N_2} a_1[k]z^{-k}} \right)
\end{align}
$$

Where individual coefficients $b_1[k],a_1[k],b_2[k],a_2[k]$ are those that represent the cascaded systems.

We want $B[k]$ and $A[k]$ that represent the cascade system, from which we can use the following Matlab/Octave function to plot the DTFT magnitude and phase responses, assumimg that the cascade is a stable LTI system so that it will have a frequency response function:

figure,freqz(B,A);

Now it can be shown by polynomial manipulations that the coefficients $B[k]$ and$A[k]$ are related to the individual coefficients $b_1[k],a_1[k],b_2[k],a_2[k]$ as the following:
$$ B[k] = b_1[k] \star b_2[k]$$ and
$$ A[k] = a_1[k] \star a_2[k]$$

Hence the required frequency response plot can be obtained in the combined call as follows:

figure,freqz(conv(b1,b2),conv(a1,a2));

where b1,a1 and b2,a2 are the coefficient vectors that represent the individual systems. The method can be generalised to N systems in cascade.


Understanding the constraint on distance between microphones in microphone array for Direction of Arrival estimation The distance between microphones in a two-microphone array for direction of arrival estimation, should be at least less than half of the wavelength of the signal that is being detected by the microphone, when the signal is periodic.

I am having trouble understanding this constraint on the distance between microphones. As far as I remember, the path difference between the sound arriving at one microphone and the other shouldn't be a multiple of half the wavelength, because if $x = \frac{\lambda}{2}$, the phase difference between the two signals for a continuous periodic signal would be zero, and thus we wouldn't be able to measure the time difference between the sound arriving at the two micropones.

However, I am unable to understand how a less than constraint appeared, and on the distance between the microphones and not on the path difference of the sound arriving at the two microphones.

Say the microphones $A$ and $B$ are omnidirectional and exactly half a wavelength apart. If the sound field propagates in the direction of the vector $A\to B$, then a chosen part of the waveform will arrive at $...,$ $A,$ $B,$ $A,$ $B,$ $...$ at intervals that equal half the period of the waveform. If the sound field propagates in the direction of the vector $B\to A$, the chosen part of the waveform will arrive at $...,$ $B,$ $A,$ $B,$ $A,$ $...$ at identical intervals. You won't be able to tell the difference between the two cases from the arrival times because for a periodic waveform there is no "first arrival". This ambiguity is resolved by the constraint you mention, as you then know that you can safely choose for comparison arrival times at $A$ and $B$ that are separated by at most some known time less than half the period of the waveform.


High modulation index PSK - carrier recovery How to recover the carrier of a PSK signal that has high modulation index (eg. 750 symbols per second at carrier frequency 1500 Hz)?  I guess that a PLL will not be able to lock onto it with such frequent phase transitions. Maybe upconverting or frequency multiplication of the signal before the PLL is the solution?

A PLL on its own will not work on the direct PSK modulation, assuming the symbols are equi-probable as that results in a nulled carrier, so there is nothing the PLL can track!  Costas-Loops are effective for BPSK and QPSK implementations, and as Dilip had suggested, for BPSK you can square your signal and then use a PLL to lock onto the 2F frequency that would come out of the squarer. Why? Because the modulation is going back and for at 0 and 180°, and squaring the signal is a frequency doubler: ($cos(\omega t)^2 = cos(0) +cos(2\omega t)$).  When you double the frequency, you also double the phase, so 0° and 180° maps to 0° and 360°, which is the same as 0° and 0° and hence you will have an unmodulated signal to lock onto, you then divide this by two and have a recovered carrier. You can do the same with QPSK but you would need to multiply by 4.   Also you must be careful with your selection of loop BW in your carrier tracking loop: If the BW is too wide (competitive with your modulation rate) you will start to track out the modulation! My rule of thumb is to use a carrier tracking loop BW that is 1/20 to 1/100th of the symbol modulation rate. If I need to optimize noise I will start with this and then model my actual implementation since this becomes a trade space with the phase noise in the LO and sampling clock jitter and tracking dynamics desired.
Here are modern all digital carrier recovery implementations that work for all PSK and QAM waveforms, based on detecting phase variation from symbol to symbol and using that to compensate for carrier offsets until the phase change is 0 (delta phase vs delta time = frequency and when 0 we are locked to the carrier):

Below shows the "Cross Product Phase Detector" conceptually, showing how taking the complex conjugate multiplication of two vectors results in a vector with a imaginary value that is proportional to the phase of the original two vectors!

This is shown more mathematically below:

We can use this concept to create a phase detector that can act directly on the modulated signal, such as the Decision Directed Phase Detector shown in the image below.  The cross hairs are the "decisions" to be made on each sample, choosing the closest constellation point (in this example QPSK), and the decided vector is complex conjugate multiplied with the vector prior to the decision, and the result is an error signal that is proportional to the phase error!

Since V1 and V2 in the image above are complex signals with I and Q, we can redraw that same implementation as in the figure below, which is the Costas Loop for a QPSK modulation!

Just add more decision thresholds and you can demodulate higher orders of modulation such as the QAM implementation below!

Of course the "lock range" is reduced as the order gets higher:

The lock range can be enhanced by using a "Symbol Aided Phase Detector" which is a common approach for acquisition where known data symbols are transmitted for channel estimation and initial capture of carrier phase and timing. Once acquired the previous methods shown can readily continue to track the carrier with random data symbols.

Below this gets all tied together showing a Decision Directed Carrier Recovery Loop as an all digital carrier tracking implementation and provides the final I Q decisions (demodulated symbols). The complete receiver would also have channel selection, AGC, equalization, timing recovery and matched filtering prior to this:


Which is the Matlab equivalent to IDL's function bytscl? From IDL code I got

openr,1,"measures.dat"

nx=200
ny=100

data=fltarr(200,100)

readf,1,data

iimage,bytscl(data)
iimage,sobel(data)

but I can't afford IDL, even isn't near easy to buy a student licence in my contry. Anyway

In a Matlab transcript I wrote:

I=importdata('measures.dat');
nx=200; ny=100;
S=reshape(I,nx,ny);

figure(1)
image(S)

B = edge(S,'Sobel',0,'both') ;

figure(2)
image(B)

So how bytscl can be done?

Thanks.

Here is my answer, thanks for the comments, this version may be wrong but anyway the image data looks better.

function [ Y ] = bytscl( X )

m = min( X( : ) );
M = max( X( : ) );
Top = 255;

    if isfloat( X )
    Fc = ( Top + 0.9999 )*( X - m )/( M - m );
    Y = cast( Fc, 'int8' );
    else
    Ic = ( Top + 1 )*( X - m - 1 )/( M - m );
    Y = cast( Ic, 'int8' );
    end

end


Fourier transform of a rect*half triangle I have to calculate the analytic expression of  Fourier transform  $$ x(t) = t{\rm rect} ( t- \frac{1}{2} ).$$ First I made the graph of these two signals and I obtained the graph I posted.

Now I wrote the analytic form of this graph that is $$ x(t) = t $$ for t from 0 to 1 and $$ x(t)=0 $$ for the rest. Now I should calculate the Fourier transform of $$ x(t) = t $$ but this should be $X(f)$ (?). My theory it’s completely wrong because my book obtained $$ X(f) = { \frac{1}{2} {\rm sinc} (f) + \frac{1}{i2\pi f}\left[({\rm sinc} (f )- \cos (\pi f ) \right]} e^{-i \pi f }. $$

There is a Fourier transform property called as "differenciation in frequency domain" which is as follows:

If the Fourier transform of $x(t)$ is $X(\omega)$, then the Fourier transform of $tx(t)$ is as below: $$\mathcal{F}(tx(t))  = j\frac{d}{d\omega} (X(\omega))$$

Where $\omega = 2\pi f$,

I think you would be able to easily derive the answer you mentioned in your question now.

Hope that helps


MATLAB fir2 - npt and lap This might not be the right place to ask this, but I'm hoping someone can explain two of the arguments in the MATLAB fir2 function.  It is a function for designing filters using the frequency sampling method.  There are two optional arguments:

npt - Number of grid points, specified as a positive integer scalar. npt must be larger than one-half the filter order: npt > n/2.
lap - Length of region around duplicate frequency points, specified as a positive integer scalar.

I am struggling to find any meaningful literature around this so if anyone could explain these two in further detail or point me to something to read that would be greatly appreciated.

npt is the number of frequency points that are used to define the desired frequency response. It's the length of the inverse FFT that is applied to the frequency domain data. It defaults to $512$ points, but if you want to design very long filters (with many taps) then you should choose a larger number.

lap defines the width of the transition band, if there is one. There is a transition band if there are two equal frequencies specified in the vector f, where a step in the desired response occurs. The wider the transition band the smaller the approximation error in the pass bands and stop bands, so there is a trade-off between steep transitions and approximation errors in the pass bands and stop bands.


Explaining a statement regarding aliasing in a document I am reading this document and in the section "Aliasing and Anti-Aliasing", the author make the following statement:

  Consider a signal with two frequency components: $f_1 = 10\,\text{Hz}$ which is
  our desired signal and $f_2 = 20\,\text{Hz}$ which is noise. Let's say that we
  sample the signal at  $f_s = 30\,\text{Hz}$. The first frequency component,
  $f_1$, will generate the following frequency components at the output
  of the multiplier, $10\,\text{Hz}, 20\,\text{Hz}, 40\,\text{Hz}, 50\,\text{Hz}, 70\,\text{Hz}$ and so on. The
  second frequency component $f_2 = 20\,\text{Hz}$ will generate the following
  frequency components at the output of the multiplier, $20\,\text{Hz}, 10\,\text{Hz},
> 50\,\text{Hz}, 40\,\text{Hz}, 80\,\text{Hz}$ and so on.

It's probably obvious to most but could someone explain where those generated frequency components come from. Thanks.

The background is that I've started reading Richard Lyons' DSP text and I really like it so far except that now I hit the bandpass sampling section and I didn't follow it. So, I've been looking for other material. If anyone knows of some good material ( books or papers ) on bandpass sampling ( I've printed out the dsp.stackexchange material and will check that out next ), it's appreciated.

The frequencies are $\left|10\,\mathrm{Hz}+k\cdot30\,\mathrm{Hz}\right|, k\in \mathbb{Z}$ and the same for $20\,\mathrm{Hz}$.

Note that real signals always have (complex conjugate) signal components for positive and negative frequencies, so a $10\,\mathrm{Hz}$ signal has components at $\pm10\,\mathrm{Hz}$.  Consequently aliasing in a $30\,\mathrm{Hz}$ grid will also reach $\mp 20\,\mathrm{Hz}$.


Why band-limit a signal? A Band-limited signal is defined as a signal whose Fourier Transform is zero above a specified frequency. Shannon Sampling Theorem is also stated for a strictly band-limited functions,i.e,

  "For a band-limited case, 2W numbers per second are sufficient."

Where, "W" is the max. frequency content.

My confusion is that why band-limit a signal? What would be the situation if a signal is not perfectly band-limited?

My understanding says that unless we have a band-limited signal, we can't get a usable "W" and aliasing would also occur. Does it make sense as an answer?

I think I know answers intuitively, but struggling to frame it in words. Please guide me.

A discrete-time signal represents a continuous-time signal that is band-limited to at most half the sampling frequency. Sampling a continuous-time signal that contains higher frequencies than half the sampling frequency will result in aliasing: Samples of those high frequencies are identical to and cannot be discerned from samples of other frequencies below half the sampling frequency. To avoid contamination of the signal by the aliased frequencies, you would first band-limit (lowpass filter) the continuous-time signal and only then sample it.


sketching frequency domain representation I've been given this signal and asked to sketch a frequency domain representation how do i do that? how do you know which frequency components are present in the signal?

This is a classical example of beat. This means that the signal (assuming that it is periodic and extends to plus and minus infinity) contains only two frequencies with the same amplitude. The frequency of the sine wave which seems to enclose the signal will be equal to the difference between the two frequencies. The frequency of the sine wave which seems to be modulated by the enclosing wave will be the mean of the two frequencies which are present in the signal.


Transform Length and event duration Does the length of a Fourier Transform affect the ability to discern a small event in the data?

This is a tricky question to ask, let me clarify. I don't care about when the event happened and I don't care about the resolution of the event (it doesn't have to be some precise measurement) but I want to SEE the event. I realize short events/pulses will affect the results in a weird ways, but it doesn't matter, I just want to know some event/pulse happened. Will the length of an FFT affect how long this event/pulse needs to be before I become aware of it?

I assume that "small event" means "short (in time) event".  Fourier transforms are not good instruments for detecting short events.  You would be better off looking for such events in the time domain, or with a wavelet transform, perhaps, that is set up specifically to detect the events you are interested in.

That being said, if for some reason the events have to be detected in the frequency domain then I would keep the transforms short.  The reason for this is if the event is short with respect to the transform length it will look like a non-ideal impulse, which basically just creates a lot of high frequency content.  Because of the small amount of energy (due to the short time) it will probably be difficult to pick the signal out of the noise.

If the transform is short, on the other hand, then its energy is a significant part of the total transform energy.  That will make it possible to pick the signal out of the noise, or even enable detection by measuring the total energy.


Implementation of a cascaded IIR filter I am attempting to implement a cascaded IIR filter in MATLAB.  I used fdatool to design the filter and exported it as an object, Hd.

In my script, I import a log file of data I captured and store my data in a vector called serial.  I then use MATLAB's filter(Hd,serial) to filter it and plot the results.  The filter works well.

In order to check and improve my understanding of how to implement this filter in real life (before I attempt to embed it using C) I want to filter my data using my own code and compare the results with those obtained using MATLAB's filter function

My Hd object consists of the following parameters:
SOSMatrix:

0.25, 0.5, 0.25, 1, -1.9513, 0.95189
0.25, 0.5, 0.25, 1, -1.9584, 0.95907
0.25, 0.5, 0.25, 1, -1.9718, 0.97246
0.25, 0.5, 0.25, 1, -1.9896, 0.99024

Scale Values:

1.4055e13
1
1
1
1

In order to implement my filter I have stored each line of my SOS matrix under the following titles:

    Section1 = [0.25 0.5 0.25 1 -1.9513 0.95189]
    Section2 = [0.25 0.5 0.25 1 -1.9584 0.95907]
    Section3 = [0.25 0.5 0.25 1 -1.9718 0.97246]
    Section4 = [0.25 0.5 0.25 1 -1.9896 0.99024]

My attempt at implementing the filter is given below (as I said earlier, my input data is called "serial"

First I initialize the outputs from each of the 4 stages to zeros...

    output_section_1=zeros(size(serial,1),1);
    output_section_2=zeros(size(serial,1),1);
    output_section_3=zeros(size(serial,1),1);
    output_section_4=zeros(size(serial,1),1);

I multiply my input (serial) by the first scale factor.  As the rest of the values are 1, I don't bother multiplying any of the other sections.

    serial=serial*1.4055e13;

Then I run a loop for my filter...

    for i=10:size(serial,1)

% Section 1
output_section_1(i)=Section1(1)*serial(i) + Section1(2)*serial(i-1) + Section1(3)*serial(i-2) - Section1(5)*serial(i-1) - Section1(6)*serial(i-2);

% Section 2
output_section_2(i)=Section2(1)*output_section_1(i) + Section2(2)*output_section_1(i-1) + Section2(3)*output_section_1(i-2) - Section2(5)*output_section_1(i-1) - Section2(6)*output_section_1(i-2);

% Section 3
output_section_3(i)=Section3(1)*output_section_2(i) + Section3(2)*output_section_2(i-1) + Section3(3)*output_section_2(i-2) - Section3(5)*output_section_2(i-1) - Section3(6)*output_section_2(i-2);

% Section 4
output_section_4(i)=Section4(1)*output_section_3(i) + Section4(2)*output_section_3(i-1) + Section4(3)*output_section_3(i-2) - Section4(5)*output_section_3(i-1) - Section4(6)*output_section_3(i-2);
    end

The output obtained using MATLAB's filter function looks good and is nice and smooth as would be expected.

The output obtained using my own implementation looks just like the input (noisy) except that the values have changed and there is an initialization at the beginning.

I'm sure I've made a very schoolboy mistake but I cannot see where I have gone wrong.

Any advice anyone could give me would be much appreciated.

You did not implement a recursive filter, but a non-recursive filter. Note that in your implementation each output only depends on possibly delayed input values, whereas in a recursive filter, the output depends on delayed input and output values.

If the coefficients of one second-order section are given by
$[b_0, b_1, b_2, 1, a_1, a_2]$, then if $x[n]$ is the input sequence, the output sequence $y[n]$ is computed as

$$y[n]=b_0x[n]+b_1x[n-1]+b_2x[n-2]-a_1y[n-1]-a_2y[n-2]$$


Simplify equation of single pole IIR transfer function Example - Consider the causal stable IIR transfer function
  $$
H(z)=\frac{K}{1-\alpha z^{-1}}, \quad 0 < \lvert \alpha\rvert 1
$$
  where $K$ and $\alpha$ are real constants
  Its square-magnitude function is given by
  $$
\lvert H\left(e^{j\omega}\right)\rvert^2 = H(z)H\left(z^{-1}\right)\bigg\vert_{z=e^{j\omega}}=\frac{K^2}{\left(1+\alpha^2\right)-2\alpha\cos\omega}
$$

Can anyone please explain how to simplify the equation to get this form? I am not able to understand it.

$$H(z)\cdot H(z^{-1})= \frac{K}{1-\alpha \cdot e^{-j \omega}} \cdot \frac{K}{1-\alpha \cdot e^{+j \omega}}$$

$$= \frac{K^2}{1-\alpha \cdot (e^{-j \omega}+e^{+j \omega})+ \alpha^2} = \frac{K^2}{1 + \alpha^2- 2 \cdot \alpha \cdot Re\{e^{j \omega}\}} = \frac{K^2}{1 + \alpha^2- 2 \cdot \alpha \cdot \cos( \omega)} $$


Simple beat detection algorithm for microcontroller For the purpose of controlling DIY disco lighting, I'm looking at implementing beat detection in a microcontroller with limited computing power (i.e. Arduino Uno). While detecting all beats would be a plus (for one, it would allow measuring BPM, which could be used to adjust the light to the style of music), it's not absolutely needed for this application since light setting must be typically changed at a much lower frequency than the beat. So what I'm really looking at is to detect "conspicuous" events in the music (beat, resume after blank, etc.) to synchronise light change, but without 100% reliability.

The constraints are the following:

Limited computing power is available. The Uno is based on a ATmega328 8-bit 16MHz MCU. If this is too constraining, I could consider moving to a Due, which is based on a AT91SAM3X8E, 32-bit 84MHz Cortex-M3 MCU.
Noisy sound input. Ideally, the sound would be captured with a small, cheap microphone. I'm not sure yet what sampling rate can be achieved with the Uno. In any case, I expect the signal to be fairly noisy. If that proved the be required, I could possibly consider implementing a proper sound input (e.g. with the rugged audio shield) but I'd rather avoid that if possible because of the constraint of hooking the music playing devices (there may be several) to the light control electronics.
Ideally, the algorithm would have a parameter to adjust sensitivity to beat event detection. That way, the sensitivity could be automatically increased to make sure events are detected even in bad conditions (low music volume, noise from conversations, etc.)

My questions are thus:

What (class of) algorithm would fit appropriately this application?
Is there any implementation available that could be readily integrated in a Arduino project?

I did one very simple bpm calculation in matlab and maybe it can help you!

source

Steps:

Decompose the signal in six frequency bands (six level) using DWT
Extrac the Envelope from all sub-bands

To make the Envelope Calculation do you need:

Apply Full Wave Rectification
Apply Low Pass Filter
Apply Downsample
Apply Mean Removal

Autocorrelation


Guassian Derivatives with orientations I am reading the paper Selective Search for Object Recognition here. In Section 3.2, they give a similarity measure between two regions of an image based on the texture of the regions with what they refer to as "fast SIFT-like measurements". On page 4, bottom right side of the page, they write:

We take Gaussian derivatives in eight orientations using

  $\sigma = 1$ for each colour channel. For each orientation for each colour channel we extract a histogram using a bin size of $10$.

I understand that a derivative of Gaussian filter is the filter of size $n \times n$ consisting of a discrete approximation of the derivative of a bivariate gaussian function of mean $0$ with some standard deviation.

What do the authors mean by "with eight orientations"? Is this some kind of modification to the filter? Any insights appreciated.

Based on Schuyler Smith's presentation on the paper, I think they calculated signed Gaussian derivative "gradients" at 4 orientations 0°, 45°, 90°, 135°, and then got from that 8 unsigned gradients at orientations 0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°.

If you have gradients $\frac{\partial}{\partial x}$ and $\frac{\partial}{\partial y}$ along coordinates $x$ and $y$ at orientations 0° and 90° (so that they are orthogonal, that is, perpendicular), you can calculate the gradient $\frac{\partial}{\partial x_\phi}$ with respect to a coordinate $x_\phi$ at any given orientation $\phi$ by:

$$\frac{\partial}{\partial x_\phi} = \cos(\phi)\frac{\partial}{\partial x} + \sin(\phi)\frac{\partial}{\partial y}.\tag{1}$$

The same way you can get the output of a Gaussian derivative filter at any orientation by taking a $\cos(\phi)$ and $\sin(\phi)$ weighted sum of the outputs of the $x$ and $y$ coordinate orientation Gaussian derivative filters.


MATLAB $\tt pwelch$ unit conversion Assuming the signal ($x$) is in volts.

The output units of pwelch are stated to be $V^2/\text{Hz}$. I wanted to confirm this so I ran:

fs = 48000;
t = [0:1/fs:5-1/fs];
x = 0.5*cos((2*pi*1000)*t);
[pxx f] = pwelch(x, 1024, 512, 1024, fs, 'onesided');
[m,i] = max(pxx);
sqrt(m*f(i))

Which results in $1.2697$ not $0.5$. There must be something I am missing! Anyone know what it is?

@Yatekii The linked you shared helped solve the problem. I needed to take into account (1) the gain of the window used and (2) summate over a short range of frequencies. Below is the updated code:

% Create signal
    fs = 48000;t = 0:1/fs:5-1/fs;
    A = 0.5;
    x = A*cos((2*pi*1000)*t);
    % Run FFT
    [pxx, f] = pwelch(x, 1024, 512, 1024, fs, 'one-sided');

    % Spectral parameters
    df = f(2) - f(1);
    CG = 0.54;
    NG = 0.3974;

    % Extract peak value
    isig = round(1000/df)+1;
    Aout = sqrt(sum(pxx(isig-5:isig+5)*df)) * sqrt(2);


Find the coefficient of an AR-model How can I infer the coefficients for an AR(2) model given an autocorrelation plot?

What I have tried so far: I can see this a cos-wave of 0.4 Hz, and then i find the values for each step $(1, 0.9211, 0.6967, 0.3624)$. Since it is an AR(2) model I have to find to coefficients ($a$ and $b$). Then I can solve two equations with two unknowns. So $a$ is $-0.9991$ and $b$ is $1.8141$. Is that the correct approach? And is it even the correct coefficients?

Given that you know the autocorrelation function, you can use the Yule-Walker equations for a generic AR model of order $p$. As you are interested in the case $p=2$, there's no need to complicate things that much.

As stated in this paper, the equations for the case of an AR(2) are given by:

$$a=\frac{r_1(1-r_2)}{1-r_1^2}$$

$$b=\frac{r_2-r_1^2}{1-r_1^2}$$

where $r_i$ are the coefficients corresponding to the autocorrelation function.


What does orthogonality mean in the context of adjacent channels in a filter bank? First off, what does orthogonality mean in the context of adjacent channels in a filter bank?

Secondly, I am reviewing some code in which the dot product of the transfer functions of adjacent channels are being used to "test orthogonality". Is this a formal method with a name? How does it work?

The mathematical definition of orthogonality between two vectors is that their dot product is zero.  It just means that there is no correlation between the two- at least at that "phase".  It is often the case that if you shifted one of the vectors you would get strong correlation.

Infinite vectors of different frequencies are always orthogonal, so in an ideal world the output of adjacent channels in a filter band would always be orthogonal.  There are two ways, though, that the real world is not ideal.

First, time limitations can introduce non-orthogonality.  The non-orthogonality is usually trivial if both of the vectors have "many" sinusoidal cycles, but can be substantial if the length is less than a cycle.

Second, non-ideal filters means that attenuated stop-band frequencies get into the filter output, which means that the adjacent channels do have frequencies in common, just at different amplitudes.


Use MATLAB to Restore a Signal from a Given Degraded Signal Using Tikhonov Regularization Anyone could share how to develop an application in MATLAB to restore the signal from a given degraded signal using Tikhonov regularization i.e restoring the signal $f$ via solving

$$
\min  || g - f \circ h ||^2  + \lambda ||\sigma \circ f ||^2
$$

where $\sigma = [1,-1]$. and $\circ$ denotes the circulant convolution operator. The setting of the value epsilon is dependent on the signal to noise ratio.

Given an audio file with noisy convoluted measurement,

The following routine are expected
Tikdeconv.m

Tikdeconv De-convolue a signal using Tikhonov regularisation J=deblurw(I,PSF,LAMBDA) de-convolute input signal $I$ by PSF using Tikhonov regularisation with regularisation parameter LAMBDA returning de-convoluted signal J. The assumption is that the signal $I$ was created by convolving a signal $J$ with kernel PSF and possible by adding /guassian white noise. You may not call or use any code from the following built-in function in MATLAB:deconvwr,deconv,deconv2

Any feedback on how to solve this would be appreciated. Thank you

The idea is to represent all operation sing Matrices.
Once it is done, it is easy to solve the problems as a Least Squares problems.

The way to represent Convolution Operation using a Matrix is by Toeplitz Matrix.
For 1D it is pretty straight forward to do (Just pay attention to boundary).

So let's take the simple model in the comment:

$$ g = f \circ h + n \rightarrow g = H f + n $$

Now, the LS solution is easy:

$$ \arg \min_{f} {\left \| H f - g \right \|}^{2}_{2} = { \left( {H}^{T} H \right) }^{-1} {H}^{T} g $$

The matrix form for the regularized model is:

$$ \arg \min_{f} {\left \| H f - g \right \|}^{2}_{2} + \lambda {\left \| \Sigma f \right \|}^{2}_{2} $$

Where $ H $ is the matrix Toeplitz Matrix generated from the Blur Kernel and $ \Sigma $ is the Toeplitz matrix generated form of the derivative operator (In this case, keeps the signal smooth to regulate the High Pass property of the first term).

This is actually Tikhonov Regularization and it is easily solved.

$$ \arg \min_{f} {\left \| H f - g \right \|}^{2}_{2} + \lambda {\left \| \Sigma f \right \|}^{2}_{2}  = { \left( {H}^{T} H + \lambda {\Sigma}^{T} \Sigma \right) }^{-1} {H}^{T} g $$

All the above is pretty straight forward to calculate in MATLAB once the the $ H $ and $ \Sigma $ matrices are formed (Given the Convolution Kernel, the function convmtx can assist you generate the matrix).

The way to create those matrices is pretty simple.
Given a signal of the length $ N $ (Samples) implies those matrices are $ N \times N $ (Assuming the output has the same number of samples).
Where the $ i $ row of the matrix applies the linear combination for the $ i $ output sample.

For instance, given the Gradient Kernel and a signal of length 5 the matrix becomes:

$$ \begin{bmatrix}
1 & -1 & 0 & 0 & 0 \\
0 & 1 & -1 & 0 & 0 \\
0 & 0 & 1 & -1 & 0 \\
0 & 0 & 0 & 1 & -1 \\
0 & 0 & 0 & 0 & 1
\end{bmatrix} $$

One must pay attention to the rows on the beginning and the end to apply the required boundary conditions.
In the case above I chose zero boundary conditions (Last row).

One should note, that in practice, the convolution matrix is not feasible due to memory limitations.
Hence it must be represented by a Sparse Matrix (Since it is mostly zeros) or the problem should be solved using a solver which doesn't require the matrix explicitly (Like Conjugate Gradient).

Another note would be that if Circular Convolution is after (Imitating DFT Based Convolution) the matrices should not only be Toeplitz but Circulant Matrices which means they enforce periodic boundary conditions (See the circular property in MATLAB padarray function).

Here is a MATLAB Code for the exercise.
I only used 800 samples due to memory constraints.

% Audio Deconvolution
% See http://dsp.stackexchange.com/questions/26433

%% General Parameters and Initialization

% The follwoing code must be executed before any of the following cells.
% The cells must be executed one after one in the order they are written.
clear();
close('all');

set(0, 'DefaultFigureWindowStyle', 'docked');
defaultLooseInset = get(0, 'DefaultAxesLooseInset');
set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

titleFontSize   = 14;
axisFotnSize    = 12;
stringFontSize  = 12;

thinLineWidth   = 2;
normalLineWidth = 3;
thickLineWidth  = 4;

smallSizeData   = 36;
mediumSizeData  = 48;
bigSizeData     = 60;

randomNumberStream = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber = 57162;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

addpath(genpath('RawData'));
addpath(genpath('AuxiliaryFunctions'));

%% Loading Data

load('input.mat');
load('handel.mat');

vInputSignal    = input;
vRefSignal      = handel;

vBlurKernel     = [0.0545; 0.2442; 0.4026; 0.2442; 0.0545];
vGradientKernel = [-1; 1]; % Forward Finite Differences
% vGradientKernel = [1; -1]; % Backward Finite Differences

%% Parameters
% Limiting the Number of Samples due to Memory Constraints
numSamplesToProcess = 800;

numSamples = size(vInputSignal, 1);

numSamples = min([numSamplesToProcess, numSamples]);

vInputSignal    = vInputSignal(1:numSamples);
vRefSignal      = vRefSignal(1:numSamples);

% Creating the Operators in Matrix Form
mBlurKernel     = convmtx(vBlurKernel, numSamples);
mGradientKernel = convmtx(vGradientKernel, numSamples);

% Enforcing size to match Convolution using 'same' property
mBlurKernel     = mBlurKernel(1:numSamples, :);
mGradientKernel = mGradientKernel(1:numSamples, :);

% Enforcing Circulant Matrix (Like DFT Based Convolution)
mBlurKernel(1, (end - 3):end) = vBlurKernel(1:4);
mBlurKernel(2, (end - 2):end) = vBlurKernel(1:3);
mBlurKernel(3, (end - 1):end) = vBlurKernel(1:2);
mBlurKernel(4, (end - 0):end) = vBlurKernel(1:1);

mGradientKernel(1, end:end) = vGradientKernel(1);

regFactor = 0.02; % Regularization Factor - Lambda

%% Least Squares Solution

vLsRestoredSignal = ((mBlurKernel.' * mBlurKernel) \ mBlurKernel.') * vInputSignal;

%% Regularized (Tikhonov) Least Squares Solution

vRegLsRestoredSignal = (((mBlurKernel.' * mBlurKernel) + (regFactor * (mGradientKernel.' * mGradientKernel)))\ mBlurKernel.') * vInputSignal;

%% Display Results

hFigure = figure();
hAxes   = axes();
hLineSeries = plot(1:numSamples, [vRefSignal, vInputSignal]);
set(hLineSeries, 'LineWidth', normalLineWidth);
set(get(hAxes, 'Title'), 'String', ['Reference and Input Signal'], ...
    'FontSize', titleFontSize);
set(get(hAxes, 'XLabel'), 'String', 'Samples Index', ...
    'FontSize', axisFotnSize);
set(get(hAxes, 'YLabel'), 'String', 'Samples Value', ...
    'FontSize', axisFotnSize);
hLegend = legend({['Reference Signal'], ['Input Signal']});

hFigure = figure();
hAxes   = axes();
hLineSeries = plot(1:numSamples, [vRefSignal, vLsRestoredSignal, vRegLsRestoredSignal]);
set(hLineSeries, 'LineWidth', normalLineWidth);
set(get(hAxes, 'Title'), 'String', ['Restored Signals and the Reference Signal'], ...
    'FontSize', titleFontSize);
set(get(hAxes, 'XLabel'), 'String', 'Samples Index', ...
    'FontSize', axisFotnSize);
set(get(hAxes, 'YLabel'), 'String', 'Samples Value', ...
    'FontSize', axisFotnSize);
hLegend = legend({['Reference Signal'], ['Least Squares Restored Signal'], ['Regularized Least Squares Restored Signal']});

hFigure = figure();
hAxes   = axes();
hLineSeries = plot(1:numSamples, [vInputSignal, vLsRestoredSignal, vRegLsRestoredSignal]);
set(hLineSeries, 'LineWidth', normalLineWidth);
set(get(hAxes, 'Title'), 'String', ['Restored Signals and the Input Signal'], ...
    'FontSize', titleFontSize);
set(get(hAxes, 'XLabel'), 'String', 'Samples Index', ...
    'FontSize', axisFotnSize);
set(get(hAxes, 'YLabel'), 'String', 'Samples Value', ...
    'FontSize', axisFotnSize);
hLegend = legend({['Input Signal'], ['Least Squares Restored Signal'], ['Regularized Least Squares Restored Signal']});

%% Restore Defaults
set(0, 'DefaultFigureWindowStyle', 'normal');
set(0, 'DefaultAxesLooseInset', defaultLooseInset);

The code is available at GitHub - GitHub StackExchangeCodes.


Capturing images on fast moving conveyor belt This may or may not be an image processing problem.
I am developing a simple system that captures images of a part as they are being manufactured and sent out along the conveyor belt.

I've solved the alignment problems, and background subtraction problem that this application requires. Everything is working in my prototype setup where I manually place parts with and without defect under the camera in different lighting conditions, position and other realistic variables. My setup uses Logitech C920 web cam, the software is able to process each frame in less than 1 second on average.

The problem arises when I start simulating the parts on the moving conveyor belt, it is unable to capture frame fast enough resulting in motion blur.

How is the issue approached in this industry?
If I switch to use a DSLR and use the correct f-stop (etc.) will that solve the blur problem? If so, wouldn't DSLR overheat since it has to capture so many frames for hours and hours and hours?

PS: speed of conveyor belt is fast but slow enough to easily track with human eyes. I don't have the information about the exact speed at this time.

Try reducing exposure time if your camera allows that. With the short exposure time the conveyor belt should be well-lit.

If that is not enough, a common solution is to use a flash which is synchronized to the camera. If a sync signal is not readily available from the camera and you can't trigger the camera with your own sync signal, you probably could find a pin on your prototype camera's circuit board that carries a signal you can sync the strobe light to.

Without synchronization, setting the strobe rate to about 2/3 of the frame rate may be a reliable solution that gives good frames separated by at most one bad frame. Whether that works depends on whether there is a rolling or a progressive shutter and on the exposure time. Other ratios closer to 1 might work better in case of a rolling shutter and/or a short exposure time.

About DSLR cameras, you should use a camera that uses an electronic shutter in the mode (still/video) you would use. Cameras do not normally overheat but a mechanical shutter is only specified for that many cycles.

Also line scanners are a possibility, but I have no idea if those are commonly used. Think of a flatbed scanner where instead of moving the scanner head the target is moving. Those never blur the image but require bright illumination along the line being imaged. That may be easier to arrange than illumination over a large area, but it's a more exotic solution than a normal area camera.

More of an anecdote: A synchronously rotating prism with parallel opposite sides, located in front of the camera, can approximately compensate for linear motion.


Difference between band pass vs. low pass + high pass? I am reading a paper in which they filter EEG signal using a notch filter to remove the power supply component of 50 Hz, then a high pass filter at 0.3 Hz to remove artifacts resulting from breathing, and finally a low pass filter with a cut-off frequency of 40 Hz.

Is there any difference to simply using a band pass filter of 0.3 – 40 Hz?

Additionally, why the notch at 50 Hz when everything 40+ Hz is thrown away, anyway?

Is there any difference to simply using a band pass filter of 0.3 – 40
  Hz?

No - but some filter design tools might not provide an option to design a band-pass filter with a unity gain "plateau" over a range of frequencies (instead, the response of the band-pass filter will be a bump at a given frequency).

  why the notch at 50 Hz when everything 40+ Hz is thrown away, anyway?

50 Hz is dangerously close to 40 Hz, and might still be in the transition band of the filter. The 40 Hz filter would have to be very steep to entirely attenuate the 50 Hz. It would be much easier to remove the 50 Hz component separately with a steep notch filter.


How we can obtain the time-integration of a sine signal? could we calculate the below integral by the Fourier series or the Fourier transform properties?

$$\int_{-\infty}^t \sin(\omega_0\tau)d\tau=?$$

No, this improper integral doesn't have a value. Since for the indefinite integral we have

$$\int\sin(\omega_0t)dt=-\frac{1}{\omega_0}\cos(\omega_0t)+C\tag{1}$$

the given definite integral would be

$$\int_{-\infty}^{t}\sin(\omega_0\tau)d\tau=-\frac{1}{\omega_0}\cos(\omega_0t)+\lim_{\tau\rightarrow -\infty}\frac{1}{\omega_0}\cos(\omega_0\tau)\tag{2}$$

but the limit in $(2)$ doesn't exist.


IIR filter relationship between number of poles, order and second order sections I have ported Iowa Hills Code Kit to a mcu to calculate IIR filter coefficients. They are used to initialize the CMSIS Biquad Cascade IIR Filters Using a Direct Form II Transposed Structure.
Trying to compare the output of my mcu program with a similar python script written using SciPy I realized that there is a discrepancy between the two, as explained below.
SciPy
Calculating a 4th order bandpass IIR filter using the following in SciPy,
sos = scipy.signal.iirfilter(4,
                            [1500, 2500],
                            btype='bandpass',
                            analog=False,
                            ftype='bessel',
                            output='sos',
                            fs=50000)

results in 8 zeros and 8 poles which is then converted to 4 Second-Order-Sections. This gives

My understanding however, is that the code above is producing an 8th order filter and that a 4th order filter should be represented as 2 second order sections and should have 4 poles and zeros. I might be missing something very obvious, but I'm not sure why signal.iirfilter produces this result!
IowaHills
To calculate a 4th order bandpass IIR fitler using IowaHills however, I would have to input numPoles=2 to get 2 second order sections as shown below,
TIIRCoeff coeffs = {0};              // filter coefficients struct (IowaHills)
TIIRFilterParams fparams;            // filter parameters struct (IowaHills)
fparams.IIRPassType = iirBPF;        // bandpass filter
fparams.ProtoType = BESSEL;          // proto bessel
fparams.BW = 0.10;                   // filter width
fparams.dBGain = 1.0;                // gain
fparams.Gamma = 0.0;                 // transition bandwith
// NumPoles = [order / 2] for BandPass and BandStop
// NumPoles = [order]     for LowPass and HighPass
fparams.NumPoles = 2;                // 4th order
fparams.OmegaC = 0.41;               // centre frequency for bandpass
fparams.Ripple = 0.0;                // only for chebyshev/elliptical
fparams.StopBanddB = 0.0;            // only for chebyshev/elliptical
coeffs = CalcIIRFilterCoeff(fparams);

this gives

Can anyone explain why the relationship between number of poles and filter order is different in the two programs? I'm assuming both of them are correct so what is the right way to interpret the outputs?

The filter order equals the number of poles, that's a fact. How people call and use parameters in software is a totally different question. Many routines that design standard IIR filters use frequency transformations to design bandpass or bandstop filters from a prototype lowpass filter. Often, the specified filter order is the order of the prototype lowpass filter, but the transformation doubles that order such that the resulting bandpass or bandstop filter has twice the order that was specified in the call of the routine. It works that way in Matlab/Octave. Maybe that's also the reason why you get twice the specified filter order with the SciPy routine and with the Iowa Hills software. Note that in both cases you obtain twice the specified filter order (or number of poles).


Is there any way of retrieving phase information with only access to real part of successive FFT frames? Due to a limitation of the WebAudio specification the FFT available in JavaScript only makes available the real part (1st half) of the calculated frequency data. Is there any way of recreating the phase information (and subsequently improving interpolation of peak frequencies) using the fluctuations in bin amplitude between frames?

webaudio isn't returning the real part of the complex FFT transform, but just the log magnitude of the FFT result.  From just this result alone, of the fluctuations between the results of multiple frames, you can't retrieve the phase information.

However, if you can get at the input data and decompose this input vector into a symmetric vector plus an anti-symmetric vector, the magnitude of the FFT of the symmetric vector will be equivalent to scaled real components of the FFT of the original input vector, and the magnitude of the FFT of the anti-symmetric vector will be equal scaled imaginary components of the FFT of the original input vector.  Thus, the decomposition plus 2 magnitude-only FFTs will allow you to determine both the magnitude and phase of the FFT components.

v_even = v + v_reversed
v_odd  = v - v_reversed


How to compute Discrete Cosine Transform Can anyone please help me understand the 2D equation of DCT?I find it hard to understand it in terms of how to apply the equation.

For example if I have an 8x8 image and I will apply 2D-DCT to every block of that image where will I start the application from left to right, top to bottom, or after left to right, top to bottom next.???

any help would be appreciated

You don't need to split image into blocks. The DCT equation can be applied to the whole image. The block division has been chosen for JPEG standard partly because DCT was costly to compute in the past (but that's not the only reason).

You can choose any size of block (including the single block, which is the image itself), then split image into the blocks and apply DCT for every block separately.

After applying DCT on, say 8x8 block (matrix), you get another 8x8 block, but with DCT coefficients (not pixel values). You normalize them, quantize them and then encode. The coefficients in the bottom right corner of the block (I think) can be stripped off to achieve compression (the article about JPEG describes this well).

It is very similar to DFT, which uses sine and cosine bases. The DCT uses cosine bases only. The DC component is located in the top-left corner of transformed image/block and in DFT the DC component is in the middle of the image/block.


How to Compensate for Frequency Offset in single carrier transmission Using Coarse and Fine Compensation? How to compensate for the frequency offset of a PAM, PSK, or QAM signal when we are using single carrier transmission, I am looking for algorithms, explanations to be able to carry out the implementation on matlab

To compensate for a frequency offset first it needs to be measured, let us take the example of  PSK.

At the trasnmitter : Suppose the signal is upconverted at the transmitter to $F_c$, then the received signal is given by $$s(t)e^{j2\pi F_ct}$$ here $s(t)$ is the baseband PSK signal.

At the Receiver: Suppose the received signal has a frequency offset, either due to doppler shift or inaccuracy of the receiver clock, let this offset be $\Delta_f$, after downconversion at the receiver, the signal sent to the ADC is $$s(t)e^{j2\pi F_ct}s^{*}(t)e^{j2\pi( F_c - \Delta_f)t}$$

$$\vert s(t) \vert^2 e^{-j2\pi\Delta_ft}$$

Where for the sake of simplicity, let us normalize $$\vert s(t) \vert^2 = 1$$

After going through the ADC the digitized samples are given by $$ y[n] = e^{j2\pi\Delta_fnT_s}$$ where $T_s$ is the sampling interval

Let us now correlate adjacent samples of this digitized signal, we get

$$ y[n+1]y[n]^* = e^{j2\pi\Delta_f(n+1)T_s}e^{-j2\pi\Delta_fnT_s}$$

$$ y[n+1]y[n]^* = e^{j2\pi\Delta_fT_s}$$

Everything other than $\Delta_f$ is known in the above equation, which can be found out using the following relation, let us denote the phase resulting from the correlation of samples as following:

$$ y[n+1]y[n]^* = u_k$$ then
 $$ \Delta_f = \frac{\Phi_k}{2\pi T_s} \tag{1}$$

Where $$\Phi_k = arctan\left( \frac{Img(u_k)}{Real(u_k)}  \right)$$

For a good result which is less impacted by noise, take average over a larger number of such correlations between incoming samples.

For the case when $\vert s(t) \vert^2$ is not unity or unit power, the value $E(s[n+1]s[n])$ is decided by the constellation. Usually a pilot is sent with a known constellation point to the receiver and trasnmitter to determine the frequency offset and channel.

For ex: in 4 QAM, you could send the constellation point/symbol $\sqrt{\frac{1}{2}} + j\sqrt{\frac{1}{2}}$

If however the SINR is low, the envelope of the signal $\vert s(t) \vert^2$, might not be predictable at the receiver in such cases it is better to work with the demodulated symbols and not the samples. The upside is better offset estimation, the downside is a longer convergence loop, because of need to demodulate and then  recover offset.

In such cases, the receiver will first try to demodulate the symbol with the frequency offset (hopefully it's not too big and can be demodulated correctly). The rotation in complex plane is directly linked to the offset as in equations above.

Now $T_s$ would be replaced by the symbol time in the above equation (1), because we would be correlating between adjacent symbols and not samples, rest of the analysis will remain exactly the same.

Offset correction at the receiver:

Once the frequency offset is estimated it can then be corrected at the receiver by applying the adjustment to the PLL at the receiver.

However, the decision to correct the offset also depends on the SINR and quality of measurements. If the SINR is low (that means the measured frequency offset is not fully trust worthy, then correct only a part of the measured frequency offset, ex : 75% in one  estimation cycle, and successful estimations will then converge to the actual offset over a period of time). So the amount of total correction once the estimation of the offset is done should be reflected in the confidence in the estimate decided by the SINR and other physical layer parameters.


BIBO stability of $1/x(t)$ I realize this is quite basic stuff but I am having trouble with the following question:

  Determine whether the following system with input $x(t)$ and output $y(t)$ is BIBO stable: $$y(t) = \frac{1}{x(t)}$$

My intuition says that it is not stable, just from thinking about the case where $x(t) = 0$. However, this would result in $y(t)$ being undefined, not going to infinity. Any help would be appreciated.

Hint: define a bound for $|x(t)|$, i.e., $|x(t)|\le A$; now try to find a positive number $B$ such that $|y(t)|\le B$ for any $|x(t)|\le A$ (that's simply the definition of BIBO stability). For the system $y(t)=1/x(t)$ it should be easy to show that the above cannot be satisfied (because $|x(t)|$ can get arbitrarily small).


Image Processing and applicability of 2D Fourier Transform As a newbie in the world of signal processing, I am having a hard time in appreciating image 2-D fourier transforms.

I am fully able to appreciate the concept of 1-D Fourier transform. Essentially, given a random causal signal, it can be decomposed into sinusoids. We cross-corelate known sinusoids (Basis functions) using FT and obtain the frequency spectrum. This is just perfect.

Please help me understand why in image processing, transformation along both axes is needed.

Given a greyscale image, like the Lena image, stored in a matrix $M \times N$, I infer the following:

Pixel intensity varies from left to right
Pixel intensity also varies top to bottom
A FT performed on 1 Row, Left to Right, gives us a spectrum of frequencies belonging to that row
Similarly, such a FT can be individually performed on each of the rows.
So in the end, we end up with M frequency lists of frequencies. e.g., Freq-List[Row-0] = {f1, f2, f4 ... fj}, Freq-List[Row-5] = {f2, f11 ..}

With this data, the Row frequency lists, will we not be able to tell how each pixel is affected by frequencies of that row? Shouldn't the row frequencies be sufficient?

Will the frequency lists along the columns also have a bearing on the pixels?

Please help me understand why in image processing, transformation along both axes is needed.

A one dimensional signal describes how does a quantity vary across, usually, time. Time, commonly represented by the symbol $t$, is the only parameter required to describe completely the signal at $t$.

A two dimensional signal describes how does a quantity vary across two parameters that are absolutely required to describe it completely. When referring to images, the quantity that is described is, usually, radiant flux. That is, how much "light" (more generally, radiation) is received by the sensor. In common handheld cameras, each pixel of an image describes how much visible light is received by some point in the scene that is viewed.

The complete set of pixels of an image, describes the variation of visible light across the surface of the camera's sensor.

When applying the Fourier Transform to a one dimensional signal, the dimension of time is transformed to a dimension of frequency and the transform breaks the signal down to a sum of sinusoids.

When applying the Fourier Transform to a two dimensional signal, its two spatial dimensions are decomposed into sums of orthogonal spatial sinusoids. To cut a long story short, if the basis functions of the Fourier Transform were not orthogonal, this trick of decomposing and recomposing would not be possible. How does this look like? It looks like a carton of eggs:

The higher the spatial frequency, the smaller the eggs (more of them fit in the same length) and vice versa.

More formally, it looks like this:

And it doesn't even have to be "symmetrical", that is, each one of its dimensions may be supporting a different spatial frequency:

In this last image, we have many more sinusoidal cycles across the $x$ dimension than across the $y$ dimension.

Therefore, whereas in the one dimensional case, a signal as complicated as the voice of a singer is decomposed into a set of simple "whistles", in two dimensions, an image as complicated as Lena is decomposed into a set of elementary little blobs. In the first case, the signal is correlated with a series of progressively increasing frequency sinusoids, in the second case, exactly the same thing happens only now, the signal is a patch of pixels and the "sinusoid" is a patch of spatial frequencies that could vary differently across the $x$ and $y$ dimension.

Now, in terms of expressing this process via the one dimensional Fourier Transform (that performs this correlation process with one set of sinusoids), the same, is applied twice.

Consider Lena. And apply the Fourier Transform across each one of its rows. What do you get? You get a set of rows in the frequency domain. These describe the set of sinusoids that visible light across the rows of the image varies. But!!! at this point, we know nothing about the set of sinusoids that describe visible light variation along the vertical.

Another way to "visualise" this, is to consider the DC bin of the row Fourier Transforms (frequency is zero). This tells you the average brightness of the pixels in each row but it still varies along the column direction! (i.e. we know nothing about the DC along the columns).

In other words, where we had an $f(x,y)$, we pass it through the Fourier Transform along the rows and we now have an $f(F_x, y)$. We are now in an intermediate state where one dimension is frequency and the other is still space.

For this reason, you apply the Fourier Transform once more along the columns of the Fourier Transform of the rows. (And to get back to the DC example, you now get one DC coefficient that describes the average brightness along the rows and the columns, that is, you get the average brightness of the image.)

Now, remember, the one dimensional Fourier Transform decomposes a signal ($x(t)$) into two series of "strength" coefficients, one for the strengths of the $\sin$ and one for the strength of the $\cos$ at different frequencies. The two dimensional Fourier Transform does exactly the same thing but now the "strength" coefficients are two dimensional. That is, some coefficient at pixel $i,j$ (in the Fourier Transformed image, i.e. the Frequency Domain) describes the contribution of a "carton of eggs" with different number of cycles in the $x$ and $y$ dimension to the whole image.

Generalising to three and higher dimensions is done similarly.

Hope this helps.

(Note, all images retrieved via Google Images and linked, rather than uploaded, to the post)


Connected Component Labeling - Equivalence Table I am trying to write a connected component labelling program, I have done the first part to create a label based on the neighbours. For simplicity I am using a 4-connectivity. I am using ANSI C and to create the equivalence table I am using two linked lists with structures defined like this:

struct T2_node
{
    size_t value;
    struct T2_node *next;
};

struct T1_node
{
     struct T2_node *T2;
     struct T1_node *next;
};

Every newly found object has a T1_node which links to a T2_node linked list that keeps the labels of the objects that particular label is connected to. As an example you can see this thresholded image (very noisy!) I am trying to find the connected components of:

A section of the equivalence table for this image is:

label id: connected labels (including its self).

835: 847, 835,
834: 846, 834,
833: 820, 833,
832: 812, 832,
831: 831,
830: 830,
829: 818, 829,
828: 828,
827: 929, 778, 840, 783, 827,
826: 826,
825: 825,
824: 993, 944, 904, 916, 839, 824,
823: 823,
822: 957, 794, 837, 822,
821: 956, 912, 888, 848, 821,
820: 1021, 997, 915, 924, 887, 812, 833, 820,
819: 819,
818: 829, 818,
817: 817,
816: 816,

A grey scale image of the labels of the same image can be seen below:

As you see, the labelled regions within a connected component do not all touch!

My final question is: What algorithm should I use to effectively merge the labels? In other words, how can I find which labels belong to one connected object? I am a graduate student in Astronomy, so I do not know the details of graph theory and set theory or how to implement them in a computer program, take C for example! (although I am interested to learn). So if there is a specific method to do this, please also introduce a reference I can learn that method from, thanks in advance.

As I've said in my comment to question, you can use function pixConnComp from leptonica library. Another variant is to use one of a common known algorithm for connected components labeling. For example you can use my realisation (from line 385) of conncomp labeling algorithm.

Most algorithms use two or three passes over image. My version have three passes:

Mark neighbours in each line: you just assign labels to each connected set of pixels line by line (labels are saving in matrix with image size).
Fill associative array with remarking: we assign a new labels to labels that we've get at step 1. For this operation we check neighbours of non-zero labels on image and assign to it a smallest value.
Rename markers due to associative array.

We can combine steps 2 and 3, but in that case algorithm would be more difficult.


Having Problem about DownConversion I hope, you are doing well. I want to ask question about baseband conversion.
Background:

I generated two AM signals in MATLAB, which one of them is shifted in time.
I multiply the signals with $\cos(2\pi f_c t)$ in order to obtain IF signal.
I filtered the IF signals with $f_c$ cut-off frequency.

In this way, I aim to downconverted to signal to baseband.
My problem is, when I try to cross-correlated the baseband signals, my result peaks at zero.
I expect to see the delay, but I could not. Am I missing something?
Notes:
I am exact sure, my cross-correlation function works perfectly. I tried this function before baseband process.
Also, in order to obtain both I and Q components, I try to multiply my signals with sin(2pifc*t). Nothing changes.
MessageSignal: Audio Signal.

If you mix up to $f_c$, but then use a filter with cut-off frequency of $f_c$, you're not downconverting your signal to baseband at all. You're just cutting off half the signal.
I think you want to re-visit the block diagram of the receiver architecture you're trying to implement. You're most definitely missing the important step, but since I'm not sure what kind of receiver you're trying to build, I can't tell you which.


What does it mean to plot a PAM signal over a number of symbols? Our course has been using the book Introduction to Communication Systems by Upamanyu Madhow for our communications course. There is a question in the book on page 87 Software lab 2.1 which is

Consider a pair of independently modulated signals, $u_c(t) = \sum_{n=1}^{n=N} b_c[n]p(t − n)$ and $u_s(t) = \sum_{n=1}^{n=N} b_s[n]p(t−n)$, where the symbols $b_c[n], b_s[n]$ are chosen with equal probability to be $+1$ and $-1$, and $p(t) = I_{[0,1]}(t)$ is a rectangular pulse. Let N = 100.

Use Matlab to plot a typical realization of $u_c(t)$ and $u_s(t)$ over 10 symbols. Make sure you sample fast enough for the plot to look reasonably “nice.”

I want to ask what is the meaning of plotting $u_c(t)$ and $u_s(t)$ over 10 symbols.

Link for the book - https://www.ece.ucsb.edu/wcsl/Publications/intro_comm_systems_madhow_jan2014b.pdf

$10$ symbols means that you choose $N=10$, i.e., you plot realizations of

$$u_c(t)=\sum_{n=0}^{9}b_c[n]p(t-n)$$

and

$$u_s(t)=\sum_{n=0}^{9}b_s[n]p(t-n)$$

in the range $t\in[0,10]$.


Why can the FFT always be mirrored in the middle of the x axis? Plotting the FFT of a signal with Matlab e.g. by the code:

y_fft = abs(fft(y));
f = fs*(0:length(y)-1)/length(y);
plot(f,y_fft);

where y is the discrete signal and 'fs' is sample rate. creates a plot where the curve can be mirrored in the middle of the x-axis e.g.

Why can the FFT always be mirrored in the middle of the x axis?

The DFT of a real signal is conjugate symmetric. For example, if your DFT result at, say, 2Hz was $1+j5$, then your DFT result at -2Hz would be $1-j5$. This is conjugate symmetry.

Of course, when you take the absolute magnitude, the result in both cases is the same, which is why you see this mirror image.

Going deeper, the reason for why a real signal can be decomposed into symmetric complex conjugate parts, is because a real signal oscillates on the real axis, and this can be thought of as the resultant of two phasors rotating around the center of the complex plane, in opposite directions.

Hope that helped.


FFTW: Logical array size in r2r transform? I'm taking a multidimensional discrete cosine transform (FFTW_REDFT11) in FFTW, and am unsure how to compute the scaling factor.  According to the documentation, taking the forward and reverse transforms results in the original image multiplied by a scaling factor:

  [C]omputing a transform followed by its inverse yields the original
  array scaled by N, where N is the logical DFT size. For REDFT00,
  N=2(n-1); for RODFT00, N=2(n+1); otherwise, N=2n.

The documentation also defines logical arrays for 1D transforms:

  [I]f you specify a size-5 REDFT00 (DCT-I) of the data abcde, it
  corresponds to the DFT of the logical even array abcdedcb of size 8. A
  size-4 REDFT10 (DCT-II) of the data abcd corresponds to the size-8
  logical DFT of the even array abcddcba, shifted by half a sample.

I'm not sure what n is in a multidimensional transform.  Say I have a 10x10 matrix.  Is n twice the total number of pixels (200)?  Or is it the number of pixels in the image if it were reflected in all dimensions (400)?

You're on the right track.  The distinction must be made between the physical size of your array, n (generally $n_0 × n_1 × n_2 × … × n_{d-1}$ except for in-place R2C transforms), and the logical size, N.  For a 10x10 matrix, the physical size is n=10*10=100, so the logical size for this particular transform is N=2n=200,  the scaling factor for the unnormalized transform-then-inverse result.

They talk a little more about that here in the Real even/odd FFTW3 docs, although they don't give much more info than the link you provided.

It's always a good idea to test with a problem that you already know the answer to,


expression for the FT of the frequency response of a system I am trying to find an expression for the Fourier Transform of the frequency response of the cascade system seen here:
Here is my approach:
$(-1)^n = (-1)^{-n}$
$v[n] = x[n]e^{j\pi n}$
$V(e^{jw}) = X(e^{j(w + \pi)})$
$w[n] = v[n] * h_1[n]$
$W(e^{jw}) = X(e^{j(w + \pi)}).H_1(e^{jw})\,\,...... (1)$
$y[n] = w[n]e^{j\pi n}$
$Y(e^{jw}) = W(e^{j(w + \pi)})\,\,...... (2)$
With equations 1 and 2, I can find $\frac{Y(e^{jw})}{X(e^{jw})}$ where I get $H(e^{jw})$ as $H_1(e^{j(w + 2\pi)})$.
Now, looking at the official solution from the book, I see:

I don't see how $Y(e^{jw}) = W(e^{j(w-\pi)})$ here. Have I made a mistake above?

Try writing your equation (2) as $Y(e^{j\omega})=W(e^{j(\omega + \pi)})=X(e^{j\omega})H_1(e^{j(\omega + \pi)})$, and now try solving for $H(e^{j\omega})$. Remember that phase shifting by $\pi$ and $-\pi$ gives the same result.


Confusion regarding discrete signal value? I know that discrete signals have finite precision. But does that means that those signals can have only integer values(for example,2,3,4 etc) but not decimal values(for example 2.1,3.2,4.2 etc)?

They can very well have decimal values, represented as either fixed point or floating point numbers.


Why in Overlap Save Method in DSP output length comes out to be same as input length? I  know that when we do convolution of $x(n)$ (Length $L$) and $h(m)$ (Length=$M$), length of output data comes out to be $L+M-1$, but when we break a large sequence of data and find output using Overlap Save Method, length comes out to be same as input.

Why do the extra $M-1$ points are appearing in this case?
In more expressive way, be $y(n) = x(n)*h(n)$:

Suppose Case 1.:

$x(n)$ (length $L$),
$h(n)$ (length $M$), leads to
$y(n)$ (length $L+M-1$)

using simple linear  convolution
Case 2.:

$x(n)$ (length $3L$),
$h(n)$ (length $M$), leads to
$y(n)$ (length $3L$)

using the Overlap Save Method that evaluats output after breaking input in 3 equal parts of
length $L$.

Your "Case 2" convolution using overlap method processing really comes out to 3L + M - 1 with M - 1 (or more) points discarded.   And another 2(M - 1) points computed and summed into the middle of those 3L points that you see.

The key concept is in the word "overlap".

So no extra points are appearing because they're always somewhere (summed or trashed).


Transform linear-phase FIR to minimum-phase FIR Given a FIR filter (non-causal) with phase zero and real coefficients given by $$H(z) = \sum_{n=-M}^{M}h[n]z^{-n}$$ with ripple $\delta_2$.

How can I obtain a filter $H_{{\rm min}}(z)$ of minimum phase of $M+1$ coefficients, from $H(z)$? How are those filter related (amplitudes, $\omega$, ripple)?

I know that if the filter is causal with real coefficients I can use this relation between the original filter and the minimum one:

$$H(z) = H_{\rm min}(z) H_{\rm ap}(z),$$

but in this exercise I can't.

What can I do?

The procedure you're looking for is called lifting and, as far as I know, it was first introduced by Hermann and Schuessler:

  O. Herrmann and H. W. Schuessler, Design of nonrecursive filters with minimum phase, Electron. Lett., 6(11): 329–330, 28th May 1970

The procedure is very well explained in this presentation by Ivan Selesnick. I'll briefly summarize it here, but it's important to look at the figures in the presentation.

Since you have a zero-phase filter, you know that $H(e^{j\omega})$ is real-valued. If you know the maximum error in the stopband $\delta_s$, you can define a new real-valued and non-negative transfer function by adding $\delta_s$ to $H(e^{j\omega})$:

$$H_2(\omega)=H(e^{j\omega})+\delta_s\ge 0\tag{1}$$

In the time domain this simply means adding $\delta_s$ to the central tap $h[0]$. In this way, all zeros on the unit circle become double zeros (look at the corresponding figure in the presentation quoted above). Now you only keep the zeros of $H_2(z)$ that are inside the unit circle, and one of each of the double zeros on the unit circle. In this way you reduce the filter order by a factor of $2$. Now you just need to scale the filter such that its passband amplitude oscillates around $1$. It is a straightforward exercise to derive the actual ripple values of the minimum phase filter given the ripple of the linear phase filter. Again, reading the presentation will help.


Variant of discrete fourier transform that isolates phase delay? I'm not quite sure of the mathematical terminology here but...

Is there a variant (or post-processing) of a Discrete Fourier Transform that separates the shape of a signal from any phase shift applied to the whole signal?

For example this image:

and this one:

should differ in at most 2 coefficients (one for each dimension in which it is shifted) in the frequency domain.

Alternatively is there a way to index images so that any image can be found in a database given a translated copy as a key?

Notes

My images are all 1D or 2D with power-of-two sizes (but not necessarily square.)
This project has nothing to do with OCR (even though I used a character glyph in the example.)  Please do not suggest OCR-specific algorithms/libraries!

Since this is a comparison of a template image, to a spatially, (and circularly) translated variant of it, then the power spectra of both will remain the same, but the phase spectra will be different.

I believe in your example, the approach of a simple circular convolution between your image, and the template will be what you want. The result is also an image, and once you have that, you search for a peak. How much off in the vertical/horizontal direction that peak is from the a peak that you would have normally got had you circularly convolved the template with itself, will give you the spatial delay you have in both horizontal and vertical space respectively.

You can implement that by taking the (same size as images) Fourier transform of your template image, (first one), and multiplying it with the Fourier transform of your image, (second one), and then inverse Fourier transforming the result back into the spatial domain. Once there, you can look for the peak as mentioned in the first paragraph.

This is assuming your images are the same size, as you stated.


Use Zero Padding to find hidden peaks I know adding zeros will not increase the frequency resolution. I have couple of waveforms in time domain measured over months. I take the fft of the waveforms to find frequency components. The problem is I get several points around each peak, but only based on visual check i found that the real peak is not the one hit by discrete frequency points.

To shovel this, 1. I can add zeros to get finer structure so that the frequency points are more close to the real peak points.

My question is: I can perform the same procedure for all waveforms. But can I really compare the peaks identified with such zero padding so that I can get the conclusion about whether there is peak shift over time or not in terms of frequency.

Zero-padding can be used to more accurately estimate the frequency of a single sinusoid in zero noise.  Since a zero-padded FFT produces a smoother looking plot, and thus shallower peaks, if there are any subtle changes in the inflections and thus peak location caused by noise, those shifts will be magnified by the zero-padded FFT, perhaps from a fraction of a FFT result bin to several bins of random shift.

So, one needs to have a good estimate of the S/N ratio to know whether and how much zero-padded will produce a usefully more accurate frequency peak estimate of a sinusoidal (or narrow band) signal.


Why is state-space representation more numerically stable than zeros-poles representation? Matlab documentation says things like:

For high order filters, the state-space form is the most numerically accurate, followed by the zero-pole-gain form. The transfer function coefficient form is the least accurate; numerical problems can arise for filter orders as low as 15.

http://www.mathworks.com/help/signal/ref/besself.html

Conversions between the TF, ZPK, and SS representations involve numerical computations and can incur loss of accuracy when overused. Because the SS and FRD representations are best suited for numerical computations, it is good practice to convert all models to SS or FRD and only use the TF and ZPK representations for construction or display purposes.

http://www.mathworks.com/products/demos/shipping/control/GSModelTypeConversions.html?product=CT#9

You can represent numeric system components using any model type. However, Numeric LTI model types are not equally well-suited for numerical computations. In general, it is recommended that you work with state-space (ss) or frequency response data (frd) models, for the following reasons:
The accuracy of computations using high-order transfer functions (tf or zpk models) is sometimes poor, particularly for MIMO or high-order systems. Conversions to a transfer function representation can incur a loss of accuracy.

http://www.mathworks.com/help/control/ug/conversion-between-model-types.html#f3-1039600
I understand why zpk is better than tf, and know that ss can handle MIMO systems that the others can't, but is ss also numerically better than zpk?  Why?  Matlab generates filter prototypes in zpk and then converts them to ss for doing filter transformations, which seems odd since you can do the filter transformations directly in zpk without converting back and forth.

This is because polynomial root finding algorithms are terrible without any implication that people code it wrongly or else. It is an ill-posed problem whereas eigenvalue computations and hence generalized eigenvalue problem solutions are much more stable.

Any mistake in any of the pole/zero data permeates to all entries of the polynomial entries and contaminates them which in turn causes big numerical inaccuracies in the polynomial root finding.

Plus, MIMO system representations in both are simply asking for trouble since common terms are very difficult to retain and after a few manipulations you can easily create dummy roots -which could have been considered to be numerical noise- becoming actual modes of the representations.

Similarly, canonical representations suffer from the same problem in state space. You can test it via creating 2 identical A matrices and perturb the entries of the row/col that holds the entries. Then create a balanced realization of the same system and perturb the A matrix entry of that representation and you'll see that the difference is usually very big.

In your comment you mention the nonuniqueness of the SS representation and it is true. But the dynamics are identical, hence they form an equivalence class up to transformations with nonzero transformation matrix $T$. That's actually a feature not a bug.


Filtering a Morse code (CW, OOK) signal I am building a morse code decoder in Python.  I am using SciPy and Numpy for FFT, filter, etc.  It is coming along ok, but the signal is noisy.  I am extracting the dominant frequency from the FFT of the input signal (which is audio from a wav file), and sampling at a rate great enough to detect whether a signal is present at these intervals.  Would I benefit from implementing a filter around this dominant frequency before processing the signal?  If so, what sort of filter, how wide, etc?  The dominant frequency is variable, but I determine it at each invocation programmatically (fft).

The bandwidth of the message content of a CW signal in Hz is roughly at least 4 times the WPM at which the dots and dashes of each Morse Code character are being keyed.  The actual signal may take up a several multiples of that bandwidth, depending on the envelope shaping or the rise/fall time of the CW transmitter used to prevent key clicks. (Plus a few harmonics of a Morse Code dot's AM modulation envelope may be required for a typical radio operator to reliably copy dots by ear, requiring a bandwidth multiple corresponding to that harmonic number).  So any audio filter bandwidth needs to be at least that wide to pass a Morse Code signal, either by measuring the WPM of the signal, or using some maximum WPM expected.


How does a low-pass filter programmatically work? I have been working on a simple low pass filter for < 100 Hz metering in my application. But so far, I am struggling with the theory behind it all. It's cool that I got it working, but I'd really enjoy it if I knew how/why it is working.

I found the following code:

void getLPCoefficientsButterworth2Pole(const int samplerate, const double cutoff, double* const ax, double* const by)
{
    double PI = M_PI;
    double sqrt2 = sqrt(2);

    double QcRaw  = (2 * PI * cutoff) / samplerate; // Find cutoff frequency in [0..PI]
    double QcWarp = tan(QcRaw); // Warp cutoff frequency

    double gain = 1 / ( 1 + sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) );

    by[2] = ( 1 - sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) ) * gain;
    by[1] = ( 2 - 2 * 2 / ( QcWarp * QcWarp ) ) * gain;
    by[0] = 1;

    ax[0] = 1 * gain;
    ax[1] = 2 * gain;
    ax[2] = 1 * gain;
}

To calculate the Coefficients. Then, in the audio samples, I 'low pass' them this way:

        xv[2] = xv[1];
        xv[1] = xv[0];

        xv[0] = pData[j];
        yv[2] = yv[1];
        yv[1] = yv[0];

        yv[0] = (ax[0] * xv[0] + ax[1] * xv[1] + ax[2] * xv[2]
                   - by[1] * yv[0]
                   - by[2] * yv[1]);

        pData[j] = yv[0];

To get a lowpass design.

I am wondering a few things:

I receive the audio samples in a simple float* array. What is that float number? The only thing I see is a number, how is that a piece of sound?
The code is using past calculations (three of them) in the new calculation per sample. Does that mean that the first 2 data samples are not filtered correctly? (not that it would matter because it's just 2 samples, but just wondering)
Trying to learn everything, I found a couple of formulas for the Butterworth (2nd Pole) filter. How are those formulas reflected in this code? None of the formulas I found have these calculations that you can see in the 'getLPCoefficientsButterworth2Pole()' function.

You asked how a low pass filter works and mentioned that the filter uses past values of you're data.  This is a non-technical discussion of what happens in a low pass filter.

The low pass filter takes differing views (shifted in time) of your signal, scales them and adds them together.  You can imagine drawing your signal 3 times, one being current, the second being shifted by one sample time, the third being shifted by 2 sample times.

At low frequencies, all the views look very similar (shifting by a single sample barely changes where you are on the signal at any instant in time).  In this case, the three versions will add together in a constructive (or at least non destructive) way, so the signal passes through the filter.

Now moving to higher frequencies, each shifted version of the signal becomes more distinct at any given instant (sample point) and in fact may even reverse in sign.  At these higher frequencies, the three versions of your signal tend to cancel out (added destructively) so the signal becomes attenuated.

Different types of filters arrange for this kind of constructive/destructive interference to take place over appropriate frequency bands to create low-pass, band-pass or high-pass filters.


Phantom harmonics when using cosine windows why do they appear and how to avoid them? Given an $L$ order cosine window, it is possible to show that the width of the main lobe is given by:
$$\omega_w = \frac{2 \pi L}{(2N+1)}$$
Where $L$ is the order of the window, $N$ is the maximum harmonic in the signal and therefore $2N+1$ is the number of samples taken over a period. If we consider the relationship
$$\Omega = w T_s$$
then the width of the window in the frequency domain is given by
$$\Omega_w = \frac{2 \pi L}{(2N+1)T_s}=\frac{2 \pi L}{T_0}$$
Where $Ts$ is the sampling time and $T_0$ is the observation interval.
The last formula clearly shows that it is possible to reduce the width of the main lobe either by oversampling by an oversampling ratio of $L$ or by increasing the observation interval by a factor of $L$.
Let's assume the following:

Sampled signal is a pure sine wave with frequency $f = \frac{1}{T}$. Then $N=1$.
The signal is sampled according to the synchronous sampling condition ($(2N+1)T_s = T_0 = T$) so no leakage error occurs.
The Hanning window is used. The width of the main lobe should be $\Omega_w = \frac{4 \pi}{(2N+1)T_s}$ since $L = 2$.

Now, if we call $P(j\Omega)$ the spectrum of the sampled signal, $P(j\Omega) = W(j\Omega) \star S(j\Omega)$ where $W$ and $S$ are the Fourier transforms of the window and of the signal, respectively and $\star$ is the convolution operator.
What would $P(j\Omega)$ look like? In the book I use, it looks like this:

STATEMENT A:
In my understanding at least, only $c$ is real (only $c$ is a replica of the real spectrum of the signal) while $a$ and $b$ are a byproduct of the convolution with the window due to its main lobe's width. In general, the book says, it is possible to say that if we are using an $L$ order cosine window, $2L-2$ phantom (or fictitious if you will) components will appear.

STATEMENT B:
If we increase $T_0$ by $L$ then the width of the Hanning window becomes the same as the one of the rectangular window ($\Omega_w = \frac{2 \pi}{(2N+1)T_s}$), therefore the fictitious harmonic components should not appear. But isn't this the same as increasing the number of samples by a ratio $L$?

I would like to understand if the two above statements are correct, in particular if the origin of the fictitious harmonic components is indeed the convolution and how one can avoid these components.
The book that I mention is Digital Signal Processing for Measurement Systems by D'Antona and Ferrero.

The magnitude and number of the lobes/harmonics is beyond my knowledge but you always obtain them when you window a signal with a finite window. No matter what windows you choose.

With some windows like the cosine family ones you have an efficient way to compute the convolution in the frequency domain, but if you compute it in the time domain the result is the same.

Hope it helps a bit


What are the possible ways to find spots with most concentrated contrast on an image? I'm thinking of pixelize + threshold, but I need a more formalized way of finding the darkest spots on the photo, with tweakable parameters.
Something I can think of is the parts with lowest spatial frequency and least luminous pixels together, i.e. the spots with darkest pixels spread over a largest area (on average) should win but small black dots should be considered too, with some sort of adaptive weighing.
desired result; left image transformed to right image:

UPD: The required result is, in fact, a metric. Which should be a ratio of locally cumulated contrast over a known area.
Local points of negligibly small size should not count, no matter how contrast they are; bigger areas with mediocre contrast features should win [because their cumulative luminosity power is higher]. The closest thing I can think of is a 2D probability distribution function; kernelized, smoothed out, spatially-oriented.

Pass a sliding window of eg 100x100 pixels over the image stepped 1 pixel at a time in each dimension. Note the maximum and the minimum pixel value. The local contrast is now the local_max/local_min.
A more robust method might be to use a tapered window (eg Gaussian) and a sliding histogram, using eg the 10-percentile and 90-percentile to avoid salt and pepper noise corrupting results.
If (as stated in your update) you are interested in contrast of features larger than a pixel, downsample your image first.
-k


Difference between these two integrators I am an electrical engineering student but new to digital signal processing. However while searching references for my senior project I came across following two different discrete integrator blocks in a MATLAB - Simulink simulation and I can't identify what these blocks do. I am familiar with continuous time integrators but this is quite new thing for me.

The above block -  Ts(z+1)/2(z-1)  is discrete time integrator with upper and lower saturation limits defined and that's all I got to know about upon clicking on it, however I want to know what does exactly that (z+1) and (z-1) means and how is it different from one shown below (KTs/(z-1)):

To the best of my knowledge above integrator is one which is direct discrete counterpart of continuous time integrator (correct me if I'm wrong), essentially Ts means it multiplies sampling time as essential in integration and 1/z-1 is inverse z-transform of a[k] = -u[-k] , right? But even if that's right, I don't understand what that term is doing here.
If my that understanding is right, what does the integrator in first block wants to tell me?
Also, the integrator in first picture is a PI torque controller, it produces reference torque command based on speed error and so that integrator is "I" block of "PI", I think but why is it different from the one in last pic. The one in last pic just... 'integrates'...much like it's continuous time counter part.
Thanks in advance and bear with my, I'm very new to this stuff.

There are different methods to approximate integration in discrete time. The most straightforward ones are the forward and backward Euler methods, and the trapezoidal method. A discrete-time system with transfer function

$$H(z)=\frac{T}{z-1}\tag{1}$$

implements the forward Euler method.

Similarly, the backward Euler method is implemented by

$$H(z)=\frac{Tz}{z-1}\tag{1}$$

Finally, the trapezoidal method is implemented by

$$H(z)=\frac{T}{2}\frac{z+1}{z-1}\tag{1}$$

These discrete-time integrators are also described on this mathworks page.


Method for calculating musical note-based magnitude spectrum My specific question is: is there a standard way to do the kind of analysis I describe below? If not, is there a useful method that I have missed?

I'm exploring ways to produce a spectral histogram of audio signal magnitude based on musical notes. I wish to calculate a spectrum spanning 10 octaves, or 120 notes. This covers a frequency range from ~20Hz to ~21,000Hz. The sample rate will be at least 44,100Hz.

The final desired result is a numerical value per note, that I can then plot on an axis of note bins. This is essentially a logarithmic x axis, with non-linear note bins. I wish to do this in realtime, on a desktop computer, and display visually at ~10Hz update rate or better.

I initially looked at the DFT (via FFT), however this has some drawbacks - the size of the DFT needs to be large enough to provide sufficient resolution at the lowest frequencies. E.g. to get ~1Hz resolution at 44100 Hz I will need to use an FFT of size 65536, which is longer than 1/10th of a second of signal (for 10Hz update). Secondly, the frequency bins for the DFT are linearly spaced in frequency, so there needs to be a warped mapping from these linear frequency bins to the non-linear note bins. In some cases many frequency bins will map to a single note bin. I understand that the shape of the frequency bins is related to the window function, so I imagine I would use something like sinc interpolation to collect frequency bins into corresponding note bins.

I've also come across some other options:

Someone claimed to have found a way to tweak the DFT to use note frequencies as the basis functions, however I wasn't able to understand how they did this and there was no math provided to back it up. Is this possible? Can the set of basis functions for the DFT be changed to match note frequencies, and if so what are the implications for scaling, phase (not so important in my application). I'm a layman but my feeling is that if the basis functions are not orthogonal then there's going to be overlap which could result in some terms counting energy more than once.

Another suggestion involves applying the Goertzel algorithm 120 times. However based on the rule-of-thumb on Wikipedia, this would be significantly more expensive than using the FFT above.

Another idea involved applying a much smaller FFT to a single octave, with filtering and downsampling to shift each octave down. This would result in 12 smaller FFTs and 11 frequency-shift/filter operations. This may be more efficient than the large FFT above.

One final point - the system can be fine-tuned in realtime, so "A440" may not actually be exactly 440Hz and may shift slightly by up to half a semitone in either direction. Therefore I need a method that can be "tuned" accordingly, or is precise enough to correctly capture energy within +/- half-semitone of a note.

Have you looked into a Chroma/pitch-class analysis? That sounds more like what you want overall (expect perhaps for the octave-folding/collapsing).

Check out librosa's CQT-based chromogram. Code to generate plot below on github:

You can use the tuning parameter to tune the musical scale, but study the full paramter set.

tuning : float
    Deviation (in cents) from A440 tuning

Finally, checkout these examples showing "enhanced chroma", to say isolate harmonic content:

https://librosa.github.io/librosa_gallery/auto_examples/plot_chroma.html#sphx-glr-auto-examples-plot-chroma-py

Regarding performance of the CQT computation, the FFT + a kernel can be used to make the computation much faster. See Judith Brown & Miller "Super" Puckette's An efficient algorithm for the calculation of a constant Q transform.


Removing vocoder effect from audio file There are many tutorials and guides explaining how to create a "robotic" voice by using a vocoder or a mix of filters that change pitch, tempo and maybe add a slight echoed delay. However, is it possible to do that opposite? How would you take a "robotic" sounding voice and make it sound more natural? I know certain effects, like adding emphasis, is still a research problem, but is it easier to remove the "tinny" sound found in a lot of synthetic voices?

I'm specifically interested in a code solution (preferably Python), because I'd like to wrap it around a text-to-speech synthesizer (like Festival or Flite) to work in real-time. Unfortunately, I'm not sure where to begin since Googling "how to make a robot voice sound more human" just returns tutorials on doing the exact opposite.

Here's an example voice of the Flite engine using an HTS voice. Especially near the end, the voice almost breaks up into a "buzz". Can that be removed, or would that actually require adding audio data?

I seriously doubt that there is a simple way to accomplish this.  Most of the work that I have seen tends to go in the direction of improving the synthesis, or using better (and more) phoneme recordings to get better quality.  If there were a simple way, I think someone would have found it by now rather than putting so much more effort into better synthesizers.

That said, you could try using a bandpass filter on the output.  The sample you referred to sounds better (at least is doesn't buzz) if you filter out everything above 800Hz, and everything below 100Hz.  That does kind of make it sound like the person speaking has a blanket stuffed in to her mouth, though.

If you would like, I could put together some python code (though you will need numpy/scipy for it) or you could use audacity or sox to see if the results are even acceptable before implementig code.