Why does the separable filter reduce the cost of computing the operator? A separable filter in image processing can be written as product of two more simple filters. Typically a 2-dimensional convolution operation is separated into 2 onedimensional filters. This reduces the cost of computing the operator.

Why is the cost of computing lower, if I use separable filter? I can not understand, why having 2 filters instead of one will increase performance

Assume you have a $N\times M$ sized image.

If you know take what is classically used, a square filter kernel, of let's say size $L\times L$, you'd need to convolve that with the picture – which gives you $N\times M$ pixels, each needing $L^2$ multiply-accumulates. So you end up with $A_{2D}=L^2MN$ operations.

Now, if you can decompose that filter into an $L$-sized horizontal and an L-sized vertical 1D-filter, you could first do all rows – that's $M$ values per row, each needing $L$ operations, so $LMN$ for all rows – and then you'd do the same with the vertical filter, so $LNM$ for all columns – and you end up with $A_{1D}=2LMN$, and you'd only need to show that

$$\begin{align}
A_{1D} &< A_{2D}\\
\iff\\
2LMN &< L^2MN &&||:LMN, \text{ legal since $L,M,N >0$}\\
\iff\\
2 &< L
\end{align}$$

most filters are larger than 2.


MATLAB $\tt butter$ function I suspect I have a very basic misconception so I'm trying to clarify that here:

A signal should be filtered before it is sampled, correct? That is, I always sample an analogue signal at some sampling frequency after it has already been filtered.

However, MATLAB's butter function acts on sampled data. See the example code here:

[b,a] = butter(6,0.6);
dataIn = randn(1000,1);
dataOut = filter(b,a,dataIn);

But once you've sampled, the aliases are mixed up with the original signal so isn't this a very problematic way of doing things?

You're right. That's why your anti-aliasing filter needs to be an analog filter. You need to restrict the bandwidth of the analog filter before sampling.

However, you can do part of the anti-aliasing filtering in the digital domain by first using a higher sampling rate and a simple analog filter, and down-sampling after applying a digital anti-aliasing filter. This is explained in the answers to this question. But in any case you'll need an analog filter first.


Filtering a QPSK signal Envelope variation due to filtering of a QPSK signal becomes more pronounced at greater phase steps.
Why is that? Why do greater phase steps (180 vs 90 degrees for example) lead to more pronounced envelope variations?

I think a picture of the QPSK constellation would be helpful to intuitively see what is occuring. Consider two possible transitions in the QPSK constellation as shown below, from "B" to "C" (180° transition) and from "A" to "C" (90° transition):

Note the difference between a 180° transition as when we move from "B" to "C" compared to a 90° transition in moving from "A" to "C". The instantaneous amplitude will be the vector magnitude at any point along the dashed lines during the transition. If the signal and has infinite bandwidth (not realistic, but if), the transition will be instant and therefore the vector magnitude will never change (will be constant envelope). However with any filtering there will be a finite time to complete the transition, and we will observe the magnitude change according to the vector length along the path taken. I placed stars where the magnitude is minimum in each case, and it is clear how in the 180° case where we pass through the origin, the magnitude approaches zero while in the 90° case the reduction in magnitude is much less.
Based on actual filtering and prior symbol history the paths will deviate somewhat from what I shown in the picture above, but the concept still holds. For example below shows a constellation of a root raised cosine signal where we can make out the same effect with the different transitions.

Minimizing the peak to average ratio is important in power sensitive applications such as Satcom and battery operated devices. This is the motivation for modified waveforms such that no trajectories go through the origin such as $\pi/4$-QPSK and waveforms that are constant envelope such as MSK and and FQPSK.


What type of circuit is responsible for convolution in the classic analog telephone? I'm interested in learning how telephones work, so I did a little bit of reading about signal processing. When I came up with the word convolution, I quickly realized the importance of this term.

To begin with, I want to know how classic analog telephones worked. The apparent simplicity of their design appeals to me.

What kind of circuit was responsible for convoluting the microphone impulses in a classic telephone receiver, so that this signal could be transmitted over wires?.

Convolution is a mathematical abstraction describing how a linear, time-invariant system affects a signal going through it.

Sometimes one explicitly designs a system to convolve a signal by a predefined impulse response (for example when building a digital filter); but more often than not, convolution is used to model various physical processes involved in a system. These physical processes can be transmission delays or dispersions, the limited bandwidth of an amplifier, transmission medium or electro-mechanical system, the explicit use of passive R, L, C network to achieve some filtering, etc... Thus, it might not be out of place to find convolution used in the description of some elements of a telephone system, but it's a modeling tool, not an actual process.

Here is an analogy: the trajectory of a cannonball is a parabola, but it makes little sense to ask which mechanical device in a cannon is responsible for computing $y = ax^2 + b$.


Harmonic Product Spectrum algorithm. I dont understand one step I was looking some Harmonic Product Spectrum algorithm examples and I came up with this:

//Implement Harmonic Product Spectrum
for(k = 0; k < buffer_size / 8; k++)
{
    sum[k] = magFFT[k] * magFFT[2*k] * magFFT[3*k];
    // find fundamental frequency (maximum value in plot)
    if( sum[k] > max_value2 && k > 0 )
    {
        max_value2 = sum[k];
        fund_freq = k;
    }
 }
 fund_freq1 = fund_freq * 8000 / buffer_size;

What I dont understand is the reason behind the last line, I tested it and it works but I couldn't find an explanation for that multiplication. Can anyone help me please?

Source of code: SOURCE

"8000 / buffer_size"  is frequency resolution, decide by  sample length or "buffer_size",

you know the sampling interval is "1/8000", the whole time length 'T' of the sequence is "buffer_size/8000", then '1/T' is frequency resolution,

in your code,if 'k' is 0, represent 0hz,is the signal's dc part,if 'k' is 1,represent '1/T' or "8000 / buffer_size" hz, then you can see 'k' correspond with frequency bin "k*8000 / buffer_size" .


Avoid saturated pixels in cross-correlation computation I am doing template maching using Normalized Cross-Correlation for measuring similarity between template and image.

However, some images from my dataset contain saturated pixels (overshoot) which should be avoided in the measurement.

How to update the cross-correlation computation to avoid saturated pixels?

I have finally found an aswer in a technical report "Image Alignment and Stitching: A Tutorial" from Richard Szeliski.

A weighted SSD (Sum of Suqared Differences) function is presented on page 13:

$$E_{WSSD}(u)=\sum_{i}w_{0}(x)w_{1}(x_{i}+u)\left[I_{1}(x_{i}+u)-I_{0}(x_{i})\right]^{2}$$

where $w_{0},w_{1}$ are spatially-varying weight functions (for template and the image). $u$ is a template offset.

The result may be divided by the overlap area to avoid bias towards smaller overlap.

This however works for SSD metric, not NCC. But it may suffice once images are normalized beforehand (using contrast/brightness adjustment computed using linear regression over corresponding pixel intensities).


Does an asymmetric Bernoulli matrix satisfy the RIP? Define an $n\times N$ sensing matrix $A$ by $A_{ij} = 0$ with probability $p$, and $A_{ij} = 1/\sqrt{n}$ with probability $1-p$.  Does $A$ satisfy the restricted isometry property?

For reference, the symmetric case is answered in the following paper:

  R.G. Baraniuk, M.A. Davenport, R.A. DeVore, and M.B. Wakin, "A simple
  proof of the restricted isometry property for random matrices,"
  Constructive Approximation, 28(3) pp. 253-263, December 2008.
  (pdf)

As others have stated in the comments, the answer is "No".  The non-zero mean of the matrix dictates that a nonzero mean vector (say, all ones), will have substantially higher gain than a random vector with zero mean (say uniformly random +1,-1).

Consider the squared norm of A times a constant vector y is expected to be n*(p*N)^2.  (iteration of expectations)

The squared norm of A times a vector x drawn uniformly from (-1,+1) is expected to be n*(p*N).  (calculable by sum of variances of Binomial distribution)

The norms of x and y are the same, but the expectation of transformed norms differ by a factor of p*N -- diverging as the dimensions grow large.

Here's matlab code to help demonstrate.

n=2000;
N=1000;
p=.9;
A=double(rand(n,N)<p);
x=sign(randn(N,1));
y=ones(N,1);
Ex_normSqAx = n*(N*p);  % E[ squared norm of A times random signs ]
Ex_normSqAy = n*(N*p)^2; % E[ squared norm of A times constant vector ]
normSqAx = norm(A*x)^2;
normSqAy = norm(A*y)^2;


Removing parallel edges from binary image This is a binary image. I want to remove the portion in the red region, and keep the one in the green region.

The red region has the property of double parallel edges, whereas green does not. If someone could suggest an algorithm it would be helpful.

Just a suggestion:

perform binary skeletonization
detect lines as connected sets of pixels
for every line, find distance to nearest pixel on some other line
compute average distance between inspected line and the nearest line
if the computed distance lays within some predefined interval (say, 3-8 pixels), consider it a double-line and remove both lines

The interval can be determined statistically because double lines have distances forming a distribution, which can be separated well from other line distances in the image.


What is the impulse response of SIMO channel and its equivalent in oversampling? 1- According to the below paragraph :

Could someone explain, What mean M here? I mean when saying with Order M ?

2- On the other hand,  When performing Fractional sampling, for example sampling with rate Ts/G (Ts is the sampling rate used at transmitter, and G the oversampling facto), we are expecting to get a discrete-time CIR with GxL channel taps at Ts/G spacing. (exactly as explained in answer below by Zeyad_Zeyad),

So the question, is it right that channel response in above equation is similar to SIMO channel response?  how ?

Thank you

Fractional sampling and SIMO are not exactly similar, SIMO system uses multi anntenna at receiver but Fractional sampling can be used with one antenna at receiver.  It's said that, Fractional sampling is often used to convert the SISO system into SIMO system without using Multi Antennas at receiver. How?

Suppose a SISO system is used, where you use single Tx antenna with single Rx antenna, Fractional sampling means you oversample the received signal by G, where G is  Oversample factor, that means your sampling rate becomes Ts/G, and Ts is the sampling rate used at transmitter. So, that means you don't change the structure of transmitter, but you oversample the receiver. How?

Suppose the received signal is :

where y,u,h are received signal, transmitted signal and channel respectively.

If y(t) is sampled with a rate of Ts/G instead of Ts, where G is the oversampling ratio, and "Ts is the sampling rate which was used at transmitter:  its polyphase components can be expressed as:

where g = 1, ..., G.

each polyphase component should be demodulated in a separated demodulator (This the one of disadvantage of using Fractional sampling). So instead of having one demodulator, we will have multi demodulators.

Coming back to the previous equation,

and the channel hg[n] is a matrix as below:

which means,

So now:

It seems similar,  but I can't confirm or give more details, I am not sure, AlexTP maybe can explain it better than me.

thank you and sorry for the notations, because I don't know how I can write the mathematical equations here.


Kalman Filter - Implementation, Parameters and Tuning First of all, this is the first time I try to make a Kalman filter.

I earlier posted the follwoing question Filter out noise and variations from speed values on StackOverflow which describes the background for this post. This is a typical sample of values I'm trying to filter. They don't necessarily have to decrease which is the case here. But the rate of change is typically like this

X------- Y
16 ---233.75
24 ---234.01
26 ---234.33
32 ---234.12
36 ---233.85
39 ---233.42
47 ---233.69
52 ---233.68
55 ---233.76
60 ---232.97
66 ---233.31
72 ---233.99

I have implemented my Kalman Filter according to this tutorial: Kalman Filter for Dummies.

My implementation looks like this (pseudocode).

//Standard deviation is 0.05. Used in calculation of Kalman gain

void updateAngle(double lastAngle){
  if(firsTimeRunning==true)
     priorEstimate = 0;               //estimate is the old one here
     priorErrorVariance = 1.2;        //errorCovariance is the old one
  else
     priorEstimate = estimate;              //estimate is the old one here
     priorErrorVariance = errorCovariance;  //errorCovariance is the old one
  rawValue = lastAngle;          //lastAngle is the newest Y-value recieved
  kalmanGain = priorErrorVariance / (priorErrVariance + 0.05);
  estimate = priorEstimate + (kalmanGain * (rawValue - priorEstimate));
  errorCovariance = (1 - kalmanGain) * priorErrVariance;
  angle = estimate;              //angle is the variable I want to update
}                                //which will be lastAngle next time

I start out with a prior estimate of 0. This seems to work fine. But what i notice is that the kalmanGain will decrease each time this update is run, which means that I trust my new values less the longer my filter is running (?).
I don't want that.

I went from just using a moving average (simple and exponential weighted) to using this. Right now I can't even get as good results as that did.

My question is if this is the right implementation and if my prior error variance and standard deviation looks good according to the sample values I have posted? My parameters are actually just picked randomly to see if I could get some good results. I have tried several different ranges but with poor results. If you have any suggestions to changes I can do, it would be really appreciated. I'm sorry if there is some obvious things missing. First time posting here too.

Kalman filters are useful when your input signal consists of noisy observations of some linear dynamical system's state. Given a series of observations of the system state, the Kalman filter aims to recursively provide better and better estimates of the underlying system's state. In order to apply it successfully, you need to have a model for the dynamics of the system whose state you're estimating. As described in detail at Wikipedia, this model describes how the state underlying system is expected to change over one time step, given its previous state, any inputs to the system, and a Gaussian-distributed stochastic component called process noise.

With that said, it's not clear from your question whether you have any such underlying model. The linked post indicated that you're consuming speed values from a sensor. These can be modeled as direct observations of a system's state (where the state is its speed), or indirect observations of its state (where the state is its position, for example). But, in order to use the Kalman framework, you need to choose a model for how that state is expected to evolve as time goes on; this additional information is used in order to generate the optimal estimate. A Kalman filter isn't a magic black box that will just "clean up" a signal that is applied to it.

With that said, the phenomenon that you alluded to, where the Kalman filter will become increasingly confident in its own output to the point where input observations become progressively ignored, does happen in practice. This can be mitigated by manually increasing the values in the process noise covariance matrix. Then, qualitatively, the model for the system's state transition contains a larger stochastic component, so the ability for the estimator to accurately predict the next state given the current state is diminished. This will reduce its reliance on its current estimate of the system state and increase its reliance upon subsequent observations, preventing the "ignoring the input" effect.


Discrete Fourier transform symmetry I was reading the chapter on discrete Fourier transforms in Lyons' book -- Understanding Digital Signal Processing -- and could not understand the last paragraph about symmetry.

  There’s an additional symmetry property of the DFT that deserves mention at this point. In practice, we’re occasionally required to determine the DFT of real input functions where the input index $n$ is defined over both positive and negative values. If that real input function is even, then $X(m)$ is always real and even; that is, if the real $x(n) = x(−n)$, then, $X_{\textrm{real}}(m)$ is in general nonzero and $X_{\textrm{imag}}(m)$ is zero. Conversely, if the real input function is odd, $x(n) = −x(−n)$, then $X_{\textrm{real}}(m)$ is always zero and $X_{\textrm{imag}}(m)$ is, in general, nonzero.

Note: $X(m) = X_{\textrm{real}}(m) + jX_{\textrm{imag}}(m)$

Firstly, what is meant by "odd" and "even"? I suspect it's the number of samples in the input signal, but that leads me to my second question,
Why is $X_{\textrm{imag}}(m)$ zero with real input functions that are even, and why, with real input functions that are odd, is $X_{\textrm{real}}(m)$ zero and $X_{\textrm{imag}}(m)$ generally non-zero?

Even & odd refer to the symmetry around $n = 0$.

Even means $x[n] = x[-n]$; you can get the part for $n < 0$ by simply mirroring the part for $n > 0$ at the $n=0$ line.

Odd means $x[n] = -x[-n]$; you can get the part for $n < 0$ by simply mirroring the part for $n > 0$ at the $n=0$ line and multiplying it by $-1$.

A cosine wave is even, sine wave is odd.

These are all just special cases of the general symmetry that says

  if it's real in one domain, it's conjugate symmetric in the other.

Conjugate symmetric means that the real part is even and the imaginary part is odd. Most people know that a real time domain signal as a conjugate symmetric spectrum, but it also goers the other way around: a conjugate symmetric time domain signal has a real valued spectrum.


These two images are very similar. Why doesn't SSIM reflect that? Here are two aligned images (greyscale, but same applies to RGB) that I want to compare using SSIM:

Clearly, these are very similar images. It seems that the first image is more blur rr However, the SSIM is only 0.52. I would expect it to be higher. Is SSIM not a good metric for measuring similarity in this case? Why is this the case?

clearly, these are very similar images

Both I and your arbitrarily picked measure for similarity would like to heartily contradict you, there.

It was you who picked SSIM, Structural SIMilarity, as measure; if that measure doesn't describe your own idea of similarity well enough, well, you might want to define what similarity is and come up with  a different measure.

Here, specifically: SSIM is especially designed to respect how focused, sharp, high-contrast images are. That's the whole point of it: although a low-pass version of your two images would indeed be very similar, the SSIM, meant as perceptive quality index, should identify that the second one is very different.

So, you picked pretty much the worst metric for your idea of similarity, and now you're "complaining" that it does what it's designed to do!


Remove high resolution spurious peaks from sinusoidal signal I have the current from a LEM transducer measured. The measurement is taken on the output of a transformer. The signal is a 50 Hz signal, measured at 100kHz. When the demand of the system is increased the current increases and therefore the amplitude of the sinusoidal gets bigger. I am looking for the maximum currents and the period over which this occurred.

Unfortunately the data is not the best quality. In some files the transducer connections appear to be bad as a signal sometimes will jump and cause a spike but the spike has many data points so using a median filter or similar is not helping me rid my system of the spikes.

I had initially applied a moving RMS filter to the signals (I have 50 different test files) and then taking the maximum of the signal to get the value and time where it occurred. However with the spikes it throws this out.

I thought about downsmapling the signal to reduce the resolution and hence the number of samples in the spikes and then doing a moving RMS but I would prefer to keep the resolution as per the original dataset.

Does anyone have any suggestions that I could try?

Since the spikes are regions of bad data you want to figure out how to identify where they are and exclude them from your analysis.  You certainly have more than enough sampling points to get a good read on the signal parameters on the regions that are good.

You indicate that you are looking for the amplitude.  What resolution is your goal, every cycle?  Are there harmonics?

I would set a DFT frame up that is two or three, maybe four cycles long.  Let's say four, for argument's sake.  That will be 8000 sample points.  Your DFT doesn't have to be nearly that dense.  A 64 point FFT should to the trick.  Just use every 8000 / 64  = 125th sample.  If you have a spikeless signal, only the bins that are multiples of four should have any significant magnitude and the bins in between will be near zero.  You can then be confident you have a clean read and find the peaks in the time domain.

If there is a spike in your DFT frame, then the in between bins will have larger values.  A little experimentation should allow you to find good threshold values.  It is important for you to frame the DFT on a whole number of cycles.  The easiest way to do this is to look for zero crossings.

Ced


Is there a relation between the dead zone of sigma delta modulator and filter order? This question refers to the link https://pure.tue.nl/ws/portalfiles/portal/1546321/521243.pdf. In chapter 4 it discusses the topic of dead zone, which appears in the delta sigma modulator. Dead zone means that when input signal to the modulator is very small than the output of the modulator will not follow the input, for small input signal the output of the modulator depends on the state of the modulator. The input must be higher than the dead zone to break the repetative pattern of the modulator.

The question: How the the dead zone is related to the filter order of the modulator? In the link it shows for first order filter but not much is written about higher order modulator and dead zone.

From the paper you link to:

  Modulators with a  higher order loop filter also exhibit a  dead-zone which can be calculated in a  similar manner. Because of the higher gain in the loop filter, the dead-zone of these modulators will be considerably smaller than the dead-zone of the first order SDM*. Very small signals (including system noise) will disturb the idle patterns of such SDMs, rendering the dead-zone undetectable in practical implementations.

And if you had a situation where system noise didn't disturb the idle patterns that cause dead zones, you could always inject some dither that's intentionally designed to have its spectral content at high frequencies.

* SDM = sigma delta modulator, in that paper's nomenclature.


Extracting Binary Magnetic-Strip Card Data from raw WAV I am faced with a tricky challenge:  To extract binary data from an iPhone magnetic strip card reader. This is what the magnetisation on the card looks like:

Source

Here is the .WAV the iPhone receives when you swipe a card (don't get your hopes up too much, it is a bonus loyalty card ;)). That's three swipes by the way, at different speeds. This is the raw SInt16 dump for the swipe I am using.

Someone seems to have done it here  but the actual data I capture isn't particularly easy to process.

The reading starts (and finishes) with an indeterminate number of 'zeros ' --  note that the wave only repeats after 2 ZEROS have been collected, this represents N-S followed by S-N:

(note that each of the three lines represents me swiping a different card; the bottom card in this image is 15 years old, so the magnetic field is clearly severely degraded in some places, not visible in this shot)

This will allow an algorithm to ascertain a clock tick.

The magnetic field reverses on each clock tick.  Also for a binary 1, the magnetic field reverses exactly in the middle of a tick:

The sequence always starts with a 1101+0(parity bit) start sentinel.   You can pick this out in all three readings in the above graph.   It is indicated more clearly in the cosmodro article I linked at the top of the question.

Here is an example of magnetic degradation (taken further along on the bottom card reading):

I am trying to figure out a sensible way to convert this waveform into its corresponding binary sequence.

I have found one PDF that goes into some detail,  but I can't figure out the algorithm they are using.

This PDF contains one interesting image:

If I could extract the red and blue lines as per this diagram,  I could use either one of them to extract the data, but I can't figure out the logic behind the construction.

So this is my question: How do I extract the binary sequence?

PS.  Note that the swipe speed is not going to be constant. So once the clock has been determined, in needs to be constantly adjusted from one tick to the next.

PPS.  Would autocorrelation catch pairs of ticks?  (seeing as ticks will alternate N-S S-N ... )

EDIT (June '12): I needed a lot of help on this one, but finally I have completed a solid reader ( http://www.magstripedecoder.com/ ).  Thanks for everyone who helped!  I recommend #musicdsp on IRC's efnet channel for anyone dedicated enough to take the challenge of getting to grips with the maths -- it's really really hard!

This is called biphase mark code, and you have to focus on the zero-crossings instead of the pulse amplitudes.  You have multiple zero crossings per pulse, though, because of the low-cut filters inherent to the pickup and the phone's mic input.  Yours drop farther than this between transitions, and cross zero:

You could restore a more pulsy shape by using a low-boost filter:

and then measure the pulse lengths by how long they spend above some threshold.  Maybe a better idea is to differentiate the input to make the transitions into large spikes, take the absolute value, and detect them when they go above some threshold:

Then measure the time between pulses, and when the time between two pulses is approximately the same as the last two pulses, it's a 0, when it's approximately half of what it was between the last two pulses, it's a 1.

The magnetic degradation you talk about should be easy to remove with a low-pass filter.


Continuous-time mathematical formula for deconvolution filters I have an elementary function $p:\mathbb{R}^2\rightarrow\mathbb{R}$ which (locally) represents an image. It's a polynomial, and its the result of the following 2D convolution:

$$p=f\star G\star \mathbf{1}_{A}$$

where $G(x,y)$ is a 2D Gaussian, and $\mathbf{1}_A$ 'box' function:

$$G(x,y)=e^{-\frac{(x^2+y^2)}{2\sigma^2}}, \quad\mathbf{1}_A(x,y)=1 \text{ for }x \text{ and } y < a, 0\text{ otherwise.}$$

Is there a closed form for a function with which I can convolve $p$, which results in $f$?

Currently I am 'avoiding' the use of such a filter by performing Fourier deconvolution, but mathematically this opens a huge can of worms.

The FT of the convolutional inverse is the point-wise inverse in the FD. So the convolutional inverse of a gaussian might have numerical problems. The FT of a gaussian is a gaussian. Its point-wise inverse thus can be large in many places.

For a box function, the situation doesn't look all that good either. Its FT is a variant on $\operatorname{sinc}$, which has many zeroes and loci of low magnitude.

You can certainly determine the convolutional inverse of a gaussian, but the inverse of a box function is not even a function. Convolving with a box will cancel any function that is periodic with a period multiple of the edge length.


Amplitude in frequency domain not as expected with discontinues time domain signal I'm having a problem where I cannot explain why the amplitude in the frequency domain not is as expected.

The scenario:

I have a time domain signal of consisting of two sine waves. on at 5Hz with amplitude 2 and one at 10Hz with amplitude 2. the signal is 10 seconds long. The 5 Hz sine wave is only present the first 5 seconds, and the 10Hz sine is only present in the last 5 seconds.

if I take an FFT of the first 5 seconds the magnitude in the frequency domain at 5Hz is as expected. And again if I take an FFT of the last 5 seconds the magnitude of the last 5 seconds is as expected.

But if I take the FFT of the entire signal both the magnitude at 5 and 10 Hz seems to have halfed in magnitude ? Is this expected or am I doing something wrong.

I'm using MathDotNet in a custom made application.

If you are seeing half the magnitude (amplitude is the height of the signal, magnitude is the size of a complex value), that means you are using a 1/N normalized DFT (BTW, my preference).  The sum will be the same (since the other half is orthogonal), you are just dividing it by twice as large of a number.  You should also find that all the even numbered bins are zero and the odd bins have values.  This is because the even bin sinusoids split evenly on your half way mark and the odd ones are integer and a half and therefore no longer orthogonal on those regions with your signal.

Hope this helps,

Ced

Followup

Demo for MBaz:

The output is:

Full: -0.000000 + i * 0.000000
Half: 230.265249 + i * 97.354586
Half Magnitude = 250.000000
Half Phase     = 0.400000

The code is:

#include <math.h>
#include <stdio.h>

//===============================================
int main( int argCount, char* argValues[] )
{

        double theSignal[1000];

        double theReal[1000];
        double theImag[1000];

        double theAngleFactor = 2.0 * M_PI / 1000;

        double theFrequency = 50.0;  // Cycles per frame

//--- Build the Arrays

        double thePhi = 0.4;

        for( int n = 0; n < 1000; n++ )
        {
            double theAngle = (double) n
                            * theAngleFactor;

            double theArg = theFrequency * theAngle;

            if( n < 500 )
            {
                theSignal[n] = cos( theArg + thePhi );
            }
            else
            {
                theSignal[n] = -cos( theArg + thePhi );
            }

            theReal[n] =  cos( theArg );
            theImag[n] = -sin( theArg );

        }

//--- Calculate the Bin Value

        double theBinReal = 0.0;
        double theBinImag = 0.0;

        for( int n = 0; n < 1000; n++ )
        {
            theBinReal += theSignal[n] * theReal[n];
            theBinImag += theSignal[n] * theImag[n];
        }

//--- Calculate for Half the Interval

        double theHalfReal = 0.0;
        double theHalfImag = 0.0;

        for( int n = 0; n < 500; n++ )
        {
            theHalfReal += theSignal[n] * theReal[n];
            theHalfImag += theSignal[n] * theImag[n];
        }

        double theMag = sqrt( theHalfReal * theHalfReal
                            + theHalfImag * theHalfImag );

        double thePhase = atan2( theHalfImag, theHalfReal );

//--- Show the Results

        printf( "Full: %f + i * %f\n" ,
                theBinReal, theBinImag );

        printf( "Half: %f + i * %f\n" ,
                theHalfReal, theHalfImag );

        printf( "Half Magnitude = %f\n", theMag );
        printf( "Half Phase     = %f\n", thePhase );
}
//===============================================

Appendum

To elaborate on my original answer.  I changed the code so the second half is defined like this, e.g. 10Hz:

                theSignal[n] = cos( 2.0 * theArg + thePhi );

Here are the sums:

Full: 230.265249 + i * 97.354586
Half: 230.265249 + i * 97.354586
Half Magnitude = 250.000000
Half Phase     = 0.400000

I left the amplitude at 1 so the magnitude of the first half unnormalized DFT is 250 (=500/2) as expected.  The sums are the same.

Epilogue - The smoking gun

Here are the unnormalized DFT calculations for the signal as specified by the OP near the bins of interest.

Bin 50 is the 5Hz bin.

Bin 100 is the 10Hz bin.

Sample rate is 100Hz.  N is 1000 for 10 seconds.

For normalized DFTs:

If you divide the left half sums by 500 (N for the left half DFT), your magnitude will be 1.0.
if you divide the whole sum by 1000 (N for the whole), your magnitude will be 0.5.

Notice that the 5Hz signal doesn't have much impact on the right side, and the 10Hz signal doesn't have much impact on the left.  As I've explained previuosly, the large values for the odd bins are due to the basis vectors in those regions being a whole integer plus a half number of cycles.

        L E F T    H A L F      R I G H T  H A L F           W H O L E
Bin       Real       Imag          Real       Imag          Real       Imag
====  =====================   =======================   =====================
 45     66.907      0.000        -7.771     -0.000        59.136      0.000
 46      0.000      0.000         0.000      0.000         0.000      0.000
 47    109.279      0.000        -7.959      0.000       101.320      0.000
 48      0.000     -0.000        -0.000     -0.000        -0.000     -0.000
 49    321.420      0.000        -8.166     -0.000       313.254      0.000
 50      0.000   -500.000         0.000      0.000         0.000   -500.000
 51   -315.264     -0.000        -8.392     -0.000      -323.656     -0.000
 52      0.000      0.000         0.000      0.000         0.000      0.000
 53   -103.118     -0.000        -8.641      0.000      -111.759     -0.000
 54      0.000     -0.000        -0.000     -0.000         0.000     -0.000
 55    -60.736      0.000        -8.915     -0.000       -69.651      0.000

        L E F T    H A L F      R I G H T  H A L F           W H O L E
Bin       Real       Imag          Real       Imag          Real       Imag
====  =====================   =======================   =====================
 95     -4.985      0.000       -65.080      0.000       -70.065      0.000
 96     -0.000     -0.000        -0.000     -0.000        -0.000     -0.000
 97     -4.714     -0.000      -107.504      0.000      -112.218      0.000
 98     -0.000      0.000         0.000     -0.000         0.000     -0.000
 99     -4.467      0.000      -319.694      0.000      -324.161      0.000
100     -0.000     -0.000         0.000   -500.000         0.000   -500.000
101     -4.240     -0.000       316.942     -0.000       312.701     -0.000
102      0.000      0.000        -0.000      0.000        -0.000      0.000
103     -4.033     -0.000       104.751     -0.000       100.718     -0.000
104     -0.000      0.000        -0.000      0.000        -0.000      0.000
105     -3.841     -0.000        62.325     -0.000        58.484     -0.000


How does both shifting and scaling a function affect the Fourier Transform? I understand how shifting and scaling work separately, but I'm unsure how to put them together.

That is to say, I understand if $$x(t) = X(\omega)$$ $$x(at) = \frac{1}{|a|}X\left(\frac{\omega}{a}\right)$$ and $$x(t-t_0) = e^{-j\omega{t_0}}X(\omega)$$ But if $x(t)$ is both scaled and shifted, as in $x(at-t_0)$, how would the Fourier Transform be determined using properties of the Fourier Transform?

What I tried was along the lines of:
$$
x(at-t_0) = x\left(a\left(t-\frac{t_0}{a}\right)\right) \\
\text{let } v(at) = x\left(a\left(t-\frac{t_0}{a}\right)\right) \\
\text{and if } v(at) \leftrightarrow \frac{1}{|a|}V\left(\frac{\omega}{a}\right) \\
\text{then } x\left(a\left(t-\frac{t_0}{a}\right)\right) = \frac{1}{|a|}V(\omega)e^{-jw\frac{t_0}{a}} \\
\therefore x(at-t_0) = \frac{1}{|a|}X\left(\frac{\omega}{a}\right)e^{-jw\frac{t_0}{a}}
$$

But I feel as though my logic rewriting $x(t)$ in terms of $v(t)$ isn't quite right and would appreciate if someone could either validate my logic or explain where my mistake is.

Your logic is correct. There are two ways to determine the answer. Yours is one of them:

$$\begin{align}v(t)=x(at)&\Longleftrightarrow \frac{1}{|a|}X(\omega/a)\\v(t-t_0/a)=x(at-t_0)&\Longleftrightarrow\frac{1}{|a|}X(\omega/a)e^{-j\omega t_0/a}\end{align}$$

The other option is

$$\begin{align}w(t)=x(t-t_0)&\Longleftrightarrow X(\omega)e^{-j\omega t_0}\\w(at)=x(at-t_0)&\Longleftrightarrow \frac{1}{|a|}X(\omega/a)e^{-j(\omega/a) t_0}\end{align}$$

Of course, both options lead to the same result.


Duality Property for DFT I was watching a youtube video for the duality property for continuous time Fourier transforms, which shows that if
Fourier transform of $x(t)$ is $X(\omega)$
then Fourier transform of $X(t)$ is $2\pi x(-\omega)$
Then how will duality look like in the case of the DFT?
Will it look like below?

If DFT of $x[n]$ is $X[k]$ then DFT of $X[n]$ is $2\pi x[-k]$

Duality in DFT would mean that if $x[n]$ has DFT coefficients as $X[k]$, then DFT of $X[n]$ would be $Nx[(N-k) \mod N]$

Proof:
Given, $$X[k] = \sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N}nk}, k=0,1,2,3,...,(N-1)$$
If we take DFT of the sequence $X[n]$, then what we get is the following :
$$Y[k] = \sum^{N-1}_{n=0}X[n]e^{-j\frac{2\pi}{N}nk} = N \left(\frac{1}{N}\sum^{N-1}_{n=0}X[n]e^{j\frac{2\pi}{N}n(-k)} \right)$$
Notice that the expression between "()" is the synthesis expression with $k^{th}$ frequency index replaced by $(N-k \mod N)^{th}$ index. Because, in DFT expression, $k$ can only take indices $0,1,2,3,...,(N-1)$, hence, we cannot have $-k$ as frequency index, but rather, $((N-k) \mod N)$
$$Y[k] = N \left(\frac{1}{N}\sum^{N-1}_{n=0}X[n]e^{j\frac{2\pi}{N}n(N-k)} \right) = Nx[(N-k) \mod N]$$
So, the way to interpret this is, you get a scaled and inverted sequence back when you take DFT of DFT, but $x[0]$ remains at $0^{th}$ index.
So, DFT of DFT of time-domain $x[n]$ gives $\{Nx[0], Nx[N-1], Nx[N-2], ..., Nx[2], Nx[1]\}$
Scaling by $N$ is the consequence of not dividing by $\frac{1}{\sqrt{N}}$ when taking DFT and incorporating this factor into IDFT expression.


Why is being jointly WSS important in signal estimation with LMMSE estimator? Assume that we have $y_k = s_k + n_k $.

We have observed $y_k$ and want to estimate $s_k$. The goal is to use LMMSE (Linear MMSE) estimator to find $y_k$ and on the other hand, we know that our filter is causal, it means we have

$$\hat{s_k} = \sum\limits_{m=-\infty}^k  h_{k-m}y_m\text,$$ but why should $s_k$ and $y_k$ be jointly WSS?

Recall that in the derivation of the linear MMSE estimator we force the error term $\hat{s}_k-s_k$ to be orthogonal to the signal $y_k$ i.e. $E[(\hat{s}_k-s_k)y_{k-j}] = 0$ for all $j$. The assumption of joint wide sense stationarity allows one to write $E[s_k y_{k-j}]$ as a function of the lag $j$, i.e. the covariance isn't a function of time $k$, just the lag. As a result, the optimal filter coefficients $h_k$ can be obtained by solving a system of linear equations.


function for sampling sine function I wanted to write a function for sampling sine function with sampling frequency sf,frequency f and length l, so I wrote the function adopting from the sampling formula, what change should I apply to my function? or is there any better way to do so?

function X = ssin(f,sf,l) %f as frequency, sf as sampling frequency, l as length

End=l/f;

t = 1/f:1/f:End;

d = 1/f:1/sf:End;

pt = pulstran(t,d,'rectpuls');

ms = sin(f.*t);

X = pt.*ms;

stem(X)

end

Is there a specific reason you need to use pulstran() ? Since you have an explicit expression for your signal, you can just sample the sine wave at the specified times as follows:

function x = ssin(f0,Fs,N)

% f0 = sine wave frequency [Hz]
% Fs = sampling frequency [Hz]
% N = number of samples

dt = 1/Fs;       % Sampling period [s]
T = N*dt;        % Record window duration [s]
t = 0:dt:(T-dt); % Time vector for sampling [s]

% generate samples at the specified times
x = sin(2*pi*f0*t); % [signal units]

stem(t,x)

end

% run the function
x = ssin(1,20,60)


Why FFT of a filtered step increases the lower values of the transform with respect to the unfiltered signal? I'm a bit new into DSP and when performing the FFT of a step filtered with a Notch filter I realized the frequency components arround the design filter frequency were attenuated (as expected), BUT the minimum frecuency increases arround tha design filter frequency. ¿WHY does this happen?
 These points are marked with an arrow and I believe it's just a matter of not fully understanding FFT. I would very much appreciate if someone could give me a hand with this.

I also include the time domain response of the filter (blue is the input, and red the output), and for more details about the filter it is defined as:

$$
N(z) = \frac{z^2 - 1.939 z + 0.9469}{z^2 - 1.868 z + 0.8764}
$$

(with a lot more decimals for better accuracy).

With a sampling period of 200 microseconds

So taking your filter and looking at its magnitude response we see that the magnitude close to zero will increase but all other magnitudes (on the other side of the dip) should be approximately unity gain (no magnification / increase in amplitude).

The fact that the DC gain of the filter is 0.94 means that the steady-state output for a step input will be 0.94, not 1.  That may be the effect that you're seeing.

Also, the notch is not very deep. You may want to redesign the filter to make it attenuate more.


How to concatenate IIR filter (Butterworth) with a Matlab function and evaluate the integral I have a function of $j\Omega$ and $v$ (velocity), i.e.:

$$A(j\Omega,v)=\frac{(j\Omega)^2}{p + \frac{j\Omega}{v}}$$

where $p$ is the pole.

I want know how $A(j\Omega,v)$ modify its behavior, when I insert another cascade filter $F(j\Omega)$.

In this case I would like study the power of this system via Parseval's theorem:

$$P(v) = \frac{1}{2 \pi} \int_{-\infty}^{+\infty} \big|A(j \Omega,v) F(j\Omega)\big|^2 \text{d} \Omega$$

Until yesterday, I have used always the same $F(j\Omega)$, so the Matlab code was long, but simple:

acc = @(w,v) ...
    abs( ...
    ((1i.*w).^2).* ...
    (1./(.1+1i.*w./v)) .*...
          ((0.0612.*(1i.*w)+1)./ ...
          (3.375e-06.*(1i.*w).^4+0.0001057.*(1i.*w).^3+0.0167.* ...
          (1i.*w).^2+0.0612.*(1i.*w)+1)) ...
    ).^2;

q = zeros(v_min, v_max);

for v = v_min:v_max;
    q(v) = integral(@(w)acc(w,v), 0, Inf)./pi;
end

Where the second part of acc is the filter.

This code works and everything is OK.

Now, I do not want put that filter, but several kind of IIR Butterworth low-pass filters.

A series, with $F_c = constant$, but with order N = 1:10.
A second series, with $N = constant$, and Fc = [1 2 4 8 16 32 64 128]Hz.

I know that the code for generate these filters is:

h  = fdesign.lowpass('N,F3dB', N, Fc, Fs);
Hd = design(h, 'butter');

but my questions is: how concatenate this filter with a Matlab function, i.e.:

A = @(w,v)(((1i.*w).^2).*(1./(.1+1i.*w./v)));

for evaluate the integral?

The functions fdesign and design are for digital filters. Use butter instead to obtain the filter coefficients of an analog butterworth filter. Then use the function freqs to evaluate the (Laplacian) transfer function of that filter at given frequencies.

To use integral I would first create a function that calculates the transfer function of the concatenation of filters A and B. For example:

function H = filter_cascade(w, v)
    % transfer function of A
    H_A = ((1i.*w).^2)./(.1+1i.*w./v);

    % the other transfer function
    [b,a] = butter(3, 0.3);
    H_B = freqs(b, a, w);

    % The transfer function of the cascade of A and B
    H = H_A .* H_B;

    % Return the magnitude square of the transfer function
    H = H.*conj(H);
end

Then pass a handle to this function to integral like this:

v = 0.1;
q = integral(@(w)filter_cascade(w,v), 0, 1)./pi;

Note that I have changed the upper limit of the integral to one.


What features describe audio signals? (Besides frequency and amplitude) I recorded sounds with a microphone and I try to distinguish them in my Java program. The frequency works quite good, but if I look at the fourier transforms it seems like there should be more features to distinguish sounds. I don´t know very much about signal processing, maybe you can help me. Here is a picture of two fourier transformations.

http://i.imgur.com/7npkt.png

I know that the frequency is determined by the index of the maximum magnitude (hope that is the right term). In the first fourier transform curve, its at 100, in the second at 12 (frequencies are 1102.5 and 132.3). But the two sounds look so differently when they are transformed, what else could I use to distinguish them?

Two remarks:

I am assuming you are plotting the real (or imaginary) part of the Fourier transform. It is much more common to work with the magnitude or squared magnitude (power spectrum).
The peak in the spectrum is a very poor measure of fundamental frequency (pitch). Take a piano note at 440 Hz, apply a notch filter to it to remove the 440 Hz component. Even if the peak in the spectrum is now at 880 Hz (roughly), it will still sound like a 440 Hz note - yes, at a frequency totally absent from the spectrum!

To answer your question, timbre is the property that distinguish sounds of same loudness and pitch. However, it is not an unidimensional characteristic with a well-defined unit. What you can extract are features such that things that "sound similar" have similar features.

Common features used for characterizing timbre include:

Spectral centroid, which indicate how "dark" or "bright" a sound will be perceived, spectral spread (a measure of bandwidth, tonality vs noisiness), and moments of order 3 and 4 (kurtosis and skewness).
Energy in each Bark band, or ratio of energies between consecutive Bark bands.
Any low-dimensionality descriptors of spectral envelope, including autoregressive coefficients, or Mel-Frequency cepstral coefficients.

See section 6 of this document.


Convolution in Spatial Domain Is Multiplication in Frequency Domain I have to prove convolution in spatial domain is equivalent to multiplication in frequency domain using two matrices.

$$
x(m, n) =
\begin{bmatrix}
1 && 2 \\
3 && 4
\end{bmatrix}
$$

$$
h(m, n) =
\begin{bmatrix}
1 && 0 \\
0 && 0
\end{bmatrix}
$$

When I used matrix method, I got the following result

$$
\begin{bmatrix}
1 && 2 && 0 \\
3 && 4 && 0 \\
0 && 0 && 0
\end{bmatrix}
$$

But when I converted both to DFT using the kernel

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

and multiply the result, I get

$$
\begin{bmatrix}
1 && 2 \\
3 && 4
\end{bmatrix}
$$

They are not equal. The dimensions are not matching. The first one is $3\times 3$ while the second one is $2\times 2$. Am I doing it correct? Is there something I am missing?

For finite size arrays, multiplication in the frequency domain is equivalent to circular convolution in the time domain, not linear convolution.  (Is that toroidal convolution in the 2D case?)

Note that with sufficient zero-padding, the results of circular convolution and linear convolution end up identical.


How do you perform OFDM modulation? I'm new to OFDM modulation. During OFDM modulation in the frequency domain, are the negative frequencies modulated with a conjugate version of the data of the positive frequencies, or do the negative and positive frequencies carry completely different data?

There are independent sub-carriers in both the positive and negative frequencies.  In other words, yes, they have completely different data.


2D FFT of Image and Sharpness I'm trying to make some simply autofocus algorithm for my device and pretty confused with 2D FFT results.
As I know for now - number of high-frequency's components in 2D FFT will say me that image is more focused, than other one where 2D FFT got more low-frequency's components. And i do not know how to calculate this parameter.

I'm thinking about applying high-pass filter, but how i should choose cut-off frequency? In 1D FFT it's pretty simple to determine just looking at the picture of 1D FFT..

What about just integrating the energy in the high frequencies?
You can do that after applying High Pass filter which is equivalent of weighting high frequencies differently.


How to calculate dB of a transmitted power why if I lose 80% power in the output of a system, Does it mean that I had lost 10×log10(0.2) = -7dB ?

If you have a system that causes 80% of its input power to be lost, then 20% of the input power will remain at the output. The power gain of the network is therefore:

$$
G = \frac{P_{out}}{P_{in}} =0.2
$$
$$
G|_{dB} = 10 \log_{10}(G) = 10 \log_{10}(0.2) \approx -7 \text{ dB}
$$


How do I create a polar plot from a multi-channel sound file? I have a multi-channel sound file from a microphone array. I want to get the polar plot in Matlab in order to get the directivity patter of the sound source. How do I do this? I did the FFT on the whole recording and I plotted the magnitude in dB but I'm afraid I'm not doing the correct thing.

Can someone tell me how I can get the polar directivity representation?
If I want to plot the directivity in terms of octave bands should I use a filter bank or can I just plot some averaged magnitudes of the frequencies that are on that specific octave band?

Sorry for the newbie questions and thanks a lot in advance.

If you consider ideal conditions in the measurements (anechoic conditions - no reflections, ideal omnidirectionality of the microphones), then you could use the recorded magnitude to represent the polar plots of your source's directivity. Most probably you would like to do it in 2D polar plots for various azimuths and elevations. You can also choose to plot either vertical or horizontal directivity patterns in this way. Alternatively, you could choose to create a "balloon plot" like the ones given by the Common Loudspeaker Format (CLF) files.

Now, let's briefly describe the process for a 2D horizontal polar plot. In the easy case where you have $N$ amount of microphones lying on the same horizontal plane, you just plot the recorded magnitude of each of them and find a way to fill in the gaps. This means that in the simplest case where they are spaced in equal angle steps around the source you could just interpolate every angle between them (by the way, if I am not mistaken, this is what happens in the CLF files, as their angular step is around 5 degrees). Keep in mind that this is NOT the most accurate/best way to do it.

Now, if you happen to either have a non-constant angular step you would have to find an appropriate formula to calculate the magnitude in between the microphone positions. In addition to that, if you happen to plot directivity patterns of horizontal planes that you either have some microphones on for some angles, but not for all of them (i.e. 3 mics for some horizontal angles lie on the plane but another $N - 3$ mics don't lie on the plane) you would also have to interpolate between adjacent microphone positions to get a magnitude value for that angle. A somewhat simple way to do it would be to use Vector Base Amplitude Panning (VBAP) (although this algorithm is used for "panning" it could very well serve as the basis for other calculations) to calculate the linear contribution of adjacent microphones to each location. Nevertheless, you would eventually have to find a way to calculate the magnitude for a $p = [x, y, z]^{T}$ point in 3D space in order to plot the magnitude on a $P = [x, y, c]^{T}$ or $P = [c, y, z]^{T}$ plane.

Now, regarding the frequency representation of your polar plots. If you consider the result of the FFT as a filter bank of linearly spaced equal bandwidth filters, then you can treat the data as a filter-banked data set and do the plotting for each frequency bin separately. This would result in most cases with a humongous amount of data and plots which wouldn't really be any useful (especially for high frequencies the human auditory system shows decreased frequency discrimination. For more info see critical bands). If you would like to have data in a "more meaningful human-related" way (this is somewhat a heavy state here, take it a bit light though) you could use filterbanks in either 1/3th-octave spaced bands (which is quite common), 1/6th-octave spaced bands (less common but a bit closer to the human auditory system pitch JND and critical bands bandwidth) or even use something like the bark scale or the mel scale to space your filters. Of course, the same applies to balloon plots.

Regarding your question about the way to "average" frequency bins to create "fractional-octave" plots, keep in mind that simple averaging, or even power/energy summation of the bins doesn't necessarily correspond to the intended "fractional octave" band of interest. For more information, you could have a look at "Generalized Fractional-Octave Smoothing of Audio and Acoustic Responses" by Hatziantoniou and Mourjopoulos.


Implementing contraharmonic filter I am at the moment trying to implement a contra harmonic filter,to remove the pepper noise from an image. But the implementation of it seems a bit incorrect, and I am not sure what is wrong.

void ContraHarmonic(Mat src, Mat dst, int kernel, double P)
{
  cout << "got it" << endl;
    Mat temp = src.clone();
    //copyMakeBorder(temp,temp,kernel-1,kernel-1,kernel-1,kernel-1,BORDER_CONSTANT,Scalar(0,0,0));

    for(int row = kernel/2; row < temp.rows - kernel/2-1; row++)
    {
        for(int col = kernel/2; col < temp.cols - kernel/2-1; col++)
        {

        double den=0,num=0;
          for(int i = -(kernel/2); i <= (kernel/2) ; i++)
          {
            for(int j = -(kernel/2) ; j <= (kernel/2); j++)
            {
                den += pow(temp.at<float>(row+i,col+j),P);
                num += pow(temp.at<float>(row+i,col+j),P+1);
            }
          }

        dst.at<float>(row,col) = num/den;
      }
    }
    dst.convertTo(dst,CV_8U);
}

The output i receive from it point at something being wrong with it.. I am not sure whether the implementation is incorrect, or the programming isn't done correctly.  This is the output i receive.. Which point at something being wrong.

--Updated--

Updated with the answer provided and the image looks like this

The problem occurred due to the pixel values of temp had the type uchar.

The solution to my problem was to output the pixel values as uchar, manipulate them as double, and store them as uchar.

double den=0,num=0;
      for(int i = -(kernel/2); i <= (kernel/2) ; i++)
      {
        for(int j = -(kernel/2) ; j <= (kernel/2); j++)
        {
            den += pow(temp.at<uchar>(row+i,col+j),P);
            num += pow(temp.at<uchar>(row+i,col+j),P+1);
        }
      }

    dst.at<uchar>(row,col) = num/den;


Algorithm for Color Balancing an Image I'm working on color-balancing images with as little input from the user of my app as possible, specifically to remove color casts caused by scanning aged paper. Photoshop's Auto Color action works great, with no input at all, but I don't even know where to begin in implementing something like that.

Where can I find an example of how to do that? I've done extensive searching and haven't found a good example. To see the correction I'm talking about, take a look at this question on Stack Overflow.

Photoshop's Auto Color is nothing more than Histogram Equalization.

It has some improvements over the classic Histogram Equalization by adding a threshold of not touching the % at the extreme values.

Try it and see for yourself:

Histogram Equalization.
Histogram Equalization (Image Processing) Part 1.
Histogram.

Good Luck.


Upsampling implementation I have problem with upsampling implementation in the following way:

My code:

def expand(x, L, t=None):
"""
This function expands a given signal L-times. It also returns new timestamps if they were given.

Parameters:
    x - Vector to expand.
    L - Expansion factor. It should be an integer.

Returns:
    y_L - Expanded signal.
    t_L - New timestamps.
"""

n =  L * len(x)

y_L = np.empty(n)
t_L = np.zeros(n)
#delta_t = np.abs(t[0]-t[1])

#delta_t_L = delta_t/L

for i in range(n):
    if i % L == 0:
        y_L[i] = x[i//L]
    else:
        y_L[i] = 0

t_L = np.arange(n)/L

#t_L = np.arange(0,len(x),delta_t_L).tolist() (Working only for [1, 2, 3, 4, 5, 6, 7, 6, 5, 4, 3, 2, 1])

if t is None:
    return y_L

return y_L, t_L

The problem is with 't_L' array. I don't know how to define it correct.

It's working for this signal '[1, 2, 3, 4, 5, 6, 7, 6, 5, 4, 3, 2, 1]' but then I tried to decimate this signal and then extend it in this way

K = 5
L = 3

y_K, t_K = decimate(x, K, t)
y_L, t_L = expand(y_K, L, t_K)

For [1, 2, 3, 4, 5, 6, 7, 6, 5, 4, 3, 2, 1], L = 2:

For [1. 6. 3.], K = 5 (Decimation factor):

Could you help me find the solution for this problem?

After upsampling with zero-insert, the higher frequency images need to be filtered out, which will serve the purpose of growing the zero insert values to the correct interpolated value in between samples. Once this is done, subsequent decimation operations can be properly done since every sample will then properly represent the signal.

For a signal that occupies BW $\in [0, \omega_1]$ with original sampling frequency at $\omega = 2\pi$ and is interpolated by $I$, after zero insert, the new sampling frequency will be at $I2\pi$ and $I-1$ images $2\omega_1$ wide will be centered at every multiple of the original sampling rate ($N2\pi$) for $N \in [1,2,\ldots I-1]$ within the new interpolated sampling rate $\omega \in [0, I2\pi]$.

In terms of the normalized angular output frequency with sampling rate $\omega = 2\pi$, the bandwidth of the original signal will be BW $\in [0, \omega_1/I]$ and the $I-1$ images will each be $2\omega_1/I$ wide and centered at multiples of $N2\pi/I$ for $N \in [1,2,\ldots I-1/I]$.

scipy.signal.firls() is a good choice for creating an interpolation filter and supports multiband filters to maximize rejection at the image frequencies. Decimation will follow a similar image folding where the $D-1$ locations that would be folded in should also be considered for filtering prior to decimation.

For further details see related post:
Choosing right cut-off frequency for a LP filter in upsampler


find phase greater than 2*pi in MATLAB I'm working on extracting the phase information in a given signal in MATLAB. I've the following vector;
signal = [exp(1i*10) exp(1i*100) exp(1i*1000) exp(1i*10000)];

When I use angle function in MATLAB to calculate the phase value, it returns the following result.
angle(signal) = -2.5664   -0.5310    0.9735   -2.8310

I want to calculate the original phase values, as follows,
  [10 100 1000 10000];

How to find the original phase values above interval [-π,π] in MATLAB?

You can, if you increase phase between samples slowly enough, using unwrap(angle(signal)). "Slowly enough" means the phase doesn't jump by more than $2 \pi$; unwrap works by tracking "total phase". Python equivalent implementation here (note you can configure discont there for jumps greater than $\pi$ (or TOL) but it's pointless with angle since it outputs within $[-\pi, \pi]$).
Otherwise, can't, as described by others. For an infinite number of inputs, you get the same output, like asking "what's $N$ in $N\ \%\ 2 = 0$" (division by 2 remainder).


Is there two Types of Scale Space in SURF algorithm? I couldn't figure out below paragraph on SURF paper and hope that someone can help me to understand it.

  Bay H., Ess A., Tuytelaars T. Van Gool L. - Speed-Up Robust Features
  (SURF), page 4, column 2

  The large scale changes, especially between the first filter within these octaves (from 9 to 15 is a change of 1:7),
  renders the sampling of scales quite crude. Therefore, we
  have also implemented a scale space with a finer sampling
  of the scales. This first doubles the size of the image, using
  linear interpolation, and then starts the first octave by filtering with a filter of size 15. Additional filter sizes are 21,
  27, 33, and 39. Then a second octave starts, again using filter which now increase their sizes by 12 pixels, after which
  a third and fourth octave follow. Now the scale change be-
  tween the first two filter is only 1.4 (21/15). The lowest
  scale for the accurate version that can be detected through
  quadratic interpolation is s = (1.2(18/9)/2)=1.2.

is that mean they used two different scale space at same time or they used each one individually and choose the one with better result?

They implemented each one separately, then compared their results to each other.

FH-9  --->  it's a detector that it's initial filter size is 9x9

FH-15 -->  it's a detector that it's initial filter size is 15x15


Relationship between amplitude and sampling rate I'm loading a signal with librosa in python. With the original sample rate of 22050 Hz, i get the following waveplot:

When i choose to resample my signal with the sample rate of 512 Hz, i get the following waveplot:

My question is, how come the amplitude is less with 512 Hz than with 22050 Hz? I assume that since the interval at which we record the amplitude is larger at 512Hz, the values should be larger too.

In order to resample without aliasing, the resampling process needs to apply an anti aliasing filter at the new Nyquist Frequency (or thereabouts). So chances are, your resampling process applied a low pass filter at 256 Hz or so and that removed A LOT of the energy.


Artificial Neural Network Preprocessing in Real-Time Applications I've come across an issue with my ANN when attempting to port my offline analysis to an online, real-time, application.

I currently train my algorithm using an array of input data, number of channels (columns) against number of samples (rows), and I preprocess the data by subtracting the samples mean and dividing through the standard deviation. This works well, and I'm getting accuracies for each of the output target values (4 targets in total) around 98%.

However, in a real-time application, where only one sample is arriving at a time, I cannot apply this same preprocessing method, as I can't determine the mean and STD on a single sample! So far this is purely a feed-forward problem, using saved weights from the offline analysis.

The obvious solution is to not use that preprocessing method, but it's the only way I can get working results. If I remove the preprocessing method all together the ANN just doesn't work, and all estimates are put in the first target value. Same goes for trying to preprocess the data based on the single sample (i.e. the same method mentioned above but determining the mean and STD on the data x axis, instead of the y axis).

I have two questions:

1) Why is preprocessing in this manner necessary? What is happening when I subtract the mean and divide through the STD that makes everything work?

2) What can I do to correct this issue? The only idea I have is to save the mean and STD value for real-time processing, and then update the mean/STD as new samples are presented.

Ideally I want a preprocessing method which I can perform on a single sample, or no preprocessing at all! Failing that, any suggestions on how you would go about this issue would be great.

Thanks.

The answer i give is probably partly help full.
What you can do is when you have a pre stored STD and MEAN(N samples long). While running the real time proces update the stored MEAN.
Create an array of you mean value. During the loop add every cycle a sample to your array and remove the oldest one. You create a kind of circular pattern. After N samples your MEAN is only the MEAN from your data.
Do the same with the STD.

Maybe you can create a different MEAN and STD from different types of data so the error is less worse

Good Luck


what exactly does scale mean in scale-space theory? In image processing, scale-space is a technique to represent an image at different scales. But this really confuses me, since I always think that the word scale just means the size of an object, so the scale-space representation of an image shall have different sizes of the same image.

However, it turns out that the scale-space representation is composed of Gaussian-convolved images with different $\sigma$, the image size doesn't change at all, why? And further, how $\sigma$ in the Gaussian function determines the scale of the convolved image ?

There's an image processing algorithm called Retinex that uses scales to perform local contrast enhancement. Here's an OK presentation on the Multi-Scale Retinex algorithm for more detail.

The image size doesn't change because you are performing convolution of a Guassian kernel with the image. In practice this is generally accomplished in the frequency domain and the image is padded to avoid the effects of circular convolution. This is standard stuff in image processing where you convolve a 'kernel' against an image. It's a just averaging surround pixels with a 2-D Gaussian function instead of something like a block of ones.

I've seen the choice of sigma be a somewhat adhoc "This looks good to me" selection process. There are ranges published in a few papers that seem to work well for most images. Essentially small values preserve the detail in the image, medium values provide a mix of general detail and global energy, and large values provide more global representations of the image energy. I think something like small is 1-25, medium 25-75, and large is +75. When you convolve and sum all of the scales, with a bit more light math, you get a local contrast enhanced image. The results are pretty good. You can try it out on some images to see the result of each to get a better feel for it.


What impact have different types of border-handling on high pass or low pass filtering? When high or low pass filters are being applied to an (grayscale) image, the borders need to be handled separately. They could be filled with 0's or 255's (b/w), repeat the nearest pixel value, or even the image could be repeated.

How do the different border-handling-methods impact the result of he filtering operations? What are the pro's and con's of using one method over an other with a HPF to improve the edges in a picture or a LPF to blur it?

OpenCV for example has various BorderTypes like BORDER_REPLICATE or BORDER_CONSTANT but I can't find any information about what happens when I use one of them with a HPF or LPF.

HPF or LPF linear operators. It means that an operation at pixel position i, j will only be affected by neighbors within the window size of your filter. Pixels outside the filter will not play a role - in other words, the effects of any filtering operation is localized around its support. Typically, if the image size post filtering is not a concern, we 0 pad an image and crop out the contents near the border with a size equal to half the filter window. To maintain the same size, it is advisable to replicate border values.


What Does the Total Variation Norm Mean in the Context of Image Processing What is the notion of total variation and how is total variation norm calculated in an image?

More exactly, I want to calculate and understand the meaning of $ \left \|X \right \|_{TV} $ if $X$ is an image.

If I optimize:
$$ min_{X} \left \|X \right \|_{TV} + {other . terms} $$ how will $X$ look like and what properties would be exaggerated in image, $X$ ?

The Total Variation of an image $ I $ can be calculated in one of 3 methods (See The Meaning of the Terms Isotropic and Anisotropic in the Total Variation Framework):

Anisotropic TV - $ \operatorname{TV} \left( I \right) = \sum_{x} \| \nabla I (x) \|_1 $.
Isotropic TV - $ \operatorname{TV} \left( I \right) = \sum_{x} \| \nabla I (x) \|_2 $.
Isotropic Squared TV - $ \operatorname{TV} \left( I \right) = \sum_{x} \| \nabla I (x) \|_2^2 $.

In practice, all yield almost the same result.
The Total Variation is a measure of how close the image is to a Piece Wise Constant function. Namely, How sparse are the changes in the gradient of the image.
You may read more at Total Variation Denoising.


Doubts and some confusion on variance for complex rv This question is in continuation of the one asked here.

Let's say that the measurement noise $w$ or any random variable is circularly Gaussian complex.

If the imaginary and real components each has a variance of 0.5, then what would be the total variance? Should I write: $w \sim CN(0,\sigma^2_w)$?
If on the other hand, if $w \sim CN(0,2\sigma^2_w)$ then does it mean that the real and imaginary components each have variance 1, so the total variance is 2. Is my understanding correct?
Is there a rule of thumb whether the variance should be 0.5 for each component or can it be anything?

If the imaginary and real components each has a variance of 0.5, then what would be the total variance?

The variance of the complex RV in that case would be $0.5 + 0.5 = 1$.

  If on the other hand, if $w∼CN(0,2σ^2_w)$ then does it mean that the real and imaginary components each has variance 1, so the total variance is 2. Is my understanding correct?

Each has variance $2σ^2_w/2 = σ^2_w$.

In other words, when you add two RV, the variance of the result is the sum of the variances of each RV.

  Is there a rule of thumb whether the variance should be 0.5 for each component?

When talking about Rayleigh fading coefficients, you want the variance per dimension to be 0.5 because that results in the average transmitted power being equal to the average received power (that is, the channel neither creates nor absorbs energy).


Male and female voice spectrum I have a set of data which consists of male and female voices. They are pronunciation of a same sentence. What's the appropriate method for getting an average spectrum for the male and female voices separately and comparing them? I can take fft of each voice but what's the next step to get an average spectrum? Also Welch's method can be applied to one voice each time. I know that this is in fact a random process and each of the voices is a realization but how can we estimate PSD using these realizations?
Edit: Thanks to Marcus Müller, here is the result which I've got using $$\frac{1}{N}( |\text{FFT}(x_1)|^2 + \dots + |\text{FFT}(x_n)|^2)$$
And here is the code:
  male = fft(k.');  %k: male voices
  male1 = (abs(male)).^2;
  male2 = sum(male1.')/36; %36: number of male voices

  female = fft(s.');  %s: female voices
  female1 = (abs(female)).^2;
  female2 = sum(female1.')/13;    %13: number of female voices

  plot(f , fftshift(male2) , f , fftshift(female2))
  legend({'male','female'},'Location','southwest')`

It would be really nice to see the other ways for estimation of PSD in this case.

in my opinion and limited experience, you can't just compare spectra of male and female voices to identify which is which.
2 decades ago i worked on two products called PurePitch and PitchDoctor that was able to shift the pitch of a voice independently of shifting the formants.  it also had the ability to increase or reduce the amount of pitch inflection of the voice (all the way to monotone).  the products had factory-supplied "presets", one was called "testosterone" and the other was called "estrogen", that could convincingly make a female voice sound masculine or a male voice sound feminine.
anyway, i don't remember the numerical parameters that were researched literally by my boss in usage of the alpha version of the product.  but it was something like 300 or 400 cents difference between the male and female voices and less than that for the difference between formants.
i think you might need a pitch detector and get lotsa words and do some statistics on the pitch to get a good machine idea of whether the pitch is or was more likely male or female.


Using superposition and time invariance to find RL circuit response I am confusing reading book by Chaparro (Signal System with Matlab) in example 2.7. He gives a question, a RL circuit with unit-step source $ v(t) = u(t) $ will result a response,

$$ i(t) = (1- e^{-t}) u(t) $$

If the input changed $ v(t) = u(t) - u(t-2) $, we are asked to the new response of RL system.

In the given solution (see the figure), he only use superposition and time invariance principle to bring the following,

$$ i(t) − i(t − 2) = 2(1 − e ^ {−t} ) u(t) − 2(1 − e ^ {(t−2)} ) u(t − 2) $$

Based on my understanding on superposition and time invariance, it should be,

$$ i(t) − i(t − 2) = (1 − e ^ {−t} )u(t) − (1 − e ^ {(-t+2)} ) u(t − 2) $$

But, I've checked it with graphical method (convolution) and Matlab script, the given answer (plot of $i(t)$, $-i(t-2)$ and $ i(t)-i(t-2)$) is true. Here is how,

First I obtain $h(t)$ from $i(t)$ by derrivation,

$$i(t) = v(t) * h(t)$$

I got $h(t)=e^{-t}u(t)$. By convolving $h(t)$ with $i(t)-i(t-2)$ I got the same plot on the book (Using matlab conv).

Anyone can explain, how superposition and time-invariant principle can bring the answer?

Thanks in advance.

Screenshot of the book (from Google Book)

Plot of $i(t)$, $-i(t-2)$ and $ i(t)-i(t-2)$

The answer to your last question is no. Nobody can explain how that answer was obtained because it's wrong. The fact that it is wrong is obvious, because due the term $e^{(t-2)}$ (with a positive power of $t$!), the output will keep increasing indefinitely, which is impossible for a stable system (such as the given RL-circuit).

Your answer is the correct one, and you must have made a mistake while checking the book's answer.


python: how to compute the sharpness features of image I am extracting the sharpness features of image as shown in the following image mentioned in a paper.

I have done with the following code.
Firstly, use the open cv convert the RGB to HSL (luminance is L  mentioned in the paper), then get L array.
and then used the Laplacian operator to get the LP. Finally, obtain the sharpness value of image based on mathematical form in the paper.

img_HLS = cv2.cvtColor(image, cv2.COLOR_BGR2HLS)
L = img_HLS[:, :, 1]
u = np.mean(L)
LP = cv2.Laplacian(L, cv2.CV_16S, ksize = 3)
s = np.sum(LP/u)

I don't know my solution is right or wrong, if there is problem, please help me correct it. Thank!

What you have implemented is not what the formula denotes, but what you have implemented will return an indication of Sharpness, possibly with some minor modifications.

The formula produces a single $\mathbb{R}$eal number as the ratio of the local sum of laplacians divided by the local mean. "Local" here refers to a square image patch (of some dimensions) centred around pixel $(x,y)$.

If you forget about the division by $\mu$ for a minute, the $LP(x,y)$ is the sum of the "local" second derivatives. This is a $\mathbb{R}$eal number.

The way you have implemented it, the sum operates over the Laplacian of the whole image.

Furthermore, as the Laplacian is a $\mathbb{R}$eal number and the image is a $\mathbb{N}$atural (bounded) number (a.k.a a non-negative integer), you may experience rounding errors. I say "may" because the matrix used by opencv is also based on integers, so the rounding comes into play at a different region.

Evaluating the Laplacian over the whole image is not so much of a "problem" if that sum you used operates across both dimensions. Usually, in platforms such as GNU Octave, MATLAB and others, a single sum over some matrix, operates across one dimension (usually it sums the columns). Therefore, your s would become a vector...not a $\mathbb{R}$eal number.

Hope this helps.


IIR Filter Design Issues : MATLAB I am designing an IIR Filter with the following designs constraints:

Passband $(0.75\pi,\pi)$
Stopband $(0,0.70\pi)$
No more than $0.5$ % ripple in the passband
No less than $40$ dB attenuation in the stopband

Using the designfilt command I am able to design the filter.

However when I am using the following code for the "Butteworth" filter I do not get any Phase response. The magnitude response is fine, but the Phase values have NaN values - "Not a Number". The same codes are able to design for "Chebyshev" and "Elliptic" filters.

Wp = 0.75;
Ws = 0.7;
Rp = 20*log10(1.005);
Rs = 40;

[Nb, Wnb] = buttord(Wp, Ws, Rp, Rs); % Nb: minimum order; Wnb: cutoff frequency
[Zb, Pb, Kb] = butter(Nb, Wnb,'high'); % Zb : zeros ; Pb : Poles , Kb : gain

SOSb = zp2sos(Zb, Pb, Kb);

freqz(SOSb)

Response for Butterworth Filter

Response for Chebyshev I filter

This would appear to be a bug in the freqz implementation. I didn't dig enough to figure out why but even this doesn't work:

freqz(SOSb(1:2,:));

Even though angle(freqz(SOSb(1:2,:))) returns a legitimate result.

But I found this as a workaround:

[h,f,s]=freqz(SOSb);
freqzplot(h,f,s);

Unfortunately freqzplot is being obsoleted so you might have to resort to:

plot(angle(freqz(SOSb)));


What is issue in square wave reconstruction? As shown in attached photo:

What is the reason that the approximate square wave $$x_N (t)$$ at right side is not like the original $$x(t)$$ at left side .?

Is the only reason is that number of coefficients N is very less(in this photo/case N=+_13 )and value of N needs to be raised ?

Well this goes to show that Fourier series is just approximation that gets more and more correct when you add more harmonics. Take a look at this:

$\dfrac{4\sin\theta}{\pi}$ is just first harmonic.

Harmonics are integer multiple of base frequency as you can see: $\sin3\theta$, $\sin5\theta$ etc. And the more of them you add to the first harmonic the more your final signal will look like a square wave.

The square wave you showed on the left will have infinitely many harmonics, and the one on the right is just a sum of first few.

GIF author Cmglee


Radar for Vital Signs – Matlab Simulation Interested to design and simulate FMCW radar for Vital Signs sensing.
The initial signal parameters are(1):

When using the following Matlab code to generate the waveform:
fc = 77e9;
c = 3e8;
tm = 50e-6;
bw = 4e9;
fs = 4e9;
waveform = phased.FMCWWaveform('SweepTime',tm,'SweepBandwidth',bw,...
    'SampleRate',fs);

sig = waveform();
subplot(211); plot(0:1/fs:tm-1/fs,real(sig));
xlabel('Time (s)'); ylabel('Amplitude (v)');
title('FMCW signal'); axis tight;
subplot(212); spectrogram(sig,32,16,32,fs,'yaxis');
title('FMCW signal spectrogram');

The results are:

When the fs is increased to 4GHz, the results are as expected.

Wondered what is the sampling frequency to be used in the Matlab simulation.
(1) https://www.diva-portal.org/smash/get/diva2:1326071/FULLTEXT02.pdf

By looking at the table, it seems to me that the sampling rate is $2 \text{ MHz}$. I came to that from $\frac{\text{Samples Per Chirp}}{\text{Chirp Duration}}$, which also matches up with what I expected that this is what the paper calls the "Fast Time Axis sampling rate". They may also show results for a $4 \text{ MHz}$ setup as well, see this excerpt from page 18 of the pdf:

A  text  file  containing  the  configurations  of  the  radar  must  be  provided to the radar at the boot up sequence. As mentioned, two different scenarios (measurement from back and front) are defined for the measurements of vital signs using the EVM. A couple of chirp parameters are different in the configurations. In the back configuration there is an ADCsampling rate of2 MHz (2000 Ksps) with 100 ADC samples per chirp, while there is an ADC sampling rate of 4 MHz (4000 Ksps) with 200 ADC samples per chirp in the front configuration. The sampling time or the chirp duration which is calculated from the division of number of ADC samples per chirp by ADC sampling rate for both configurations is 50 μs.

I read this to be the conventional sample rate that everybody thinks of. They also have the "Slow Time Axis sampling rate", which is how often they take a frame of samples and feed that into their processing. I think I'd call this something different, like "Frame Rate", since it really is independent of how many samples are taken per second. All this is explained on page 18-19.


Number of datapoints after decimation in discrete wavelet wavelet analysis I am following an example provided by MATLAB on 1-D wavelet decomposition (DWT). The number of datapoints is 10000. If we do a level 3 decomposition following a similar example provided by MATLAB using https://www.mathworks.com/help/wavelet/ref/wavedec.html

The length of the coefficients "l" after wavelet decomposition is

[127 127 252 501 1000], corresponding to the length of approximation coefficients at level 3, number of detail coefficients at level 3, number of details coefficients at level 2 and number of detail coefficients at level 1. Why are these numbers not exactly divided by 2?

For example, 1000/2 = 500, it is 501. Similarly at next level 500/2 should be 250 but it is 252 and at the third level 250/2 should be 125 but it is 127.

Secondly, every textbook says that DWT needs 2^n data points. Even if we choose a number such as 8192, the number of decimated coefficients are not exactly and successively divided by 2. What is the reason behind this?

Thanks.

As with Fourier, as with IIR filtering, discrete wavelet forward implementation, and how they can be inversed, depends of lot on how you treat the samples that you don't know, for standard finite support data. In other words: as you convolve, what do you do on the left side and the right side of the signal.

There is a huge literature on how to perform a non-expansive decomposition, that generates the same number of coefficients that the number of samples. But yo can see the constraints easily: when I have 3 samples, how to subsample this evenly into a low-pass and a high-pass series?

To made a long story short, if the signal length can be divided by $2^L$, $L$ wavelet levels can be obtained in general with the appropriate signal extension, governed by Matlab's dwtmode: "the options The DWT associated with the symmetric, smooth, zero, and periodic extension modes are slightly redundant."
This is generally due to how wavelet filters overlap where the signal is not defined. This is related to the convolution with lowpass and highpass filters as well, before decimation. And convolution tends to increase the signal's length. A filter with support $K$ convolved with a signal of length $N$ generally yields a filtered signal of length $N+K-1$. So if you decimate and when to keep all necessary information, you need to keep a little more coefficients at each level. There are nice description in Chapter 4 of "Wavelets in chemistry", Beata Walczak, 2000.

However in some case, like with the Haar wavelet, or as detailed in C. Brislawn, 1996 Classification of Nonexpansive Symmetric Extension Transforms for Multirate Filter Banks, nonexpansive schemes can be designed.

In practice, people who really want to be nonexpansive uses modified DWT algorithms, like the lifting scheme, or wavelets on the interval, that are not 100% pure discrete wavelet transforms.


Filtering in frequency and time domain I am trying to filter a sinusoid ($f_{sample}=50 Hz$) ,, i need to remove all frequencies above 10 Hz with accepted transition region 2 Hz (stop band freq=12 Hz)

the question is why after filtering this signal is the 20 Hz still appearing in the FFT? is there something wrong ? I know results should be the same in frequency and time domain filtering but the first in multiplication and the 2nd in convolution.

here is the code :

 t=[0:1/50:1];
 y=cos(2*pi*20*t)+0.5*sin(2*pi*5*t);
 fsamp=50;
 frange=[10 12];
 mag=[1 0];
 devs=[0.02 0.02];
 frange=[10 12];
 [n,fo,ao,w] = firpmord(frange,mag,devs,fsamp);
 n = n + rem(n,2);
 b=firpm(n,fo,ao,w,'hilbert');
 %%% filter in time domain
 final1=filter(b,1,y);
 finalfft=abs(fft(final1));
 t=[0:1/50:(length(finalfft)-1)/50];
 f=[(t/t(end))*50];
 figure(1);plot(f,abs(fft(y)));
 figure(2);plot(f,finalfft);

 %%% filter in freq domain
 paddedfft=fft(y);
 bfft=fft(b,length(y));
 finalfft=paddedfft.*bfft;
 final1=ifft(finalfft);
 t=[0:1/50:(length(finalfft)-1)/50];
 f=[(t/t(end))*50];
 t=[0:1/50:(length(final1)-1)/50];
 figure(3);plot(f,abs(fft(final1)));

Some time ago I wrote some codes just to test some filters methods, and i could see how can be hard make certain frequencies disappear completely, some methods cause ripple effects in the signal and when you apply FFT on the filtred signal to see the frequency response, you note that the frequency you want to eliminate still be there, attenuated, but you can still see!

Lets go, three(3)  sinusoidal concatenated at 250 hz, 400 hz and 900hz sampled in 8000hz and I will try to eliminate 900hz from the original signal!

My First test, AKA a very simple equalizer in Frequency Domain, calculate the dB attenuation 10.^(dB/20) and multiply by the FFT bin that you want eliminate, for dB attenuation -24dB my result was:

My Second test EQ biquad filter by Robert Bristow-Johnson (time domain), for dB attenuation -24dB my result was:

My Third test biquad LPF by Robert Bristow-Johnson (time domain), my result was:

And my last test using Sinc Filter (Frequency Domain), my result was:

And that's all, in short, everything depends on what you really need !


Questions related to Laplace Transform While studying Laplace transform, I also some questions which I want to understand:

a) We used to say that Laplace transform include both real and imaginary part whereas in Fourier transform we only have imaginary part. But when we have to say about convergence we also choose Real part to be either >0  or <0 . I want to know why we ignore imaginary part?

b) If we have any function x(t)  how do we determine that we have to take bilateral integral or unilateral integral. In the above case we have u(t)  with the function so our limits are changed. But if we don't have $u(t)$ function with the exponential like $$ x(t)=e^{at}$$ than how can we select bilateral or unilateral

Q2: what will be the Laplace Transform of $$f(t) = e^{at}$$

The Fourier transform is the Laplace transform along the imaginary axis in the complex plane.

The convergence of the Laplace transform ignores the complex part, as the imaginary part breaks down the signal into sinusoids: which are bounded, and so have no effect on the convergence.


Cross Correlation, how can any signals except the trivial cases be uncorrelated? I have been working on a problem where I am trying to subtract two Power Spectral Densities (PSD's) from one another in the following way. $$|F(k) - G(k)|^2 = |F(k)|^2 - |G(k)|^2$$

Which is only valid if the cross correlation terms are 0.

It's my understanding that for signals $f$ and $g$, the cross correlation is given by
$$(f \star g)(\tau)  = \int_{-\infty}^{\infty} \overline{f(t)}g(t+\tau) dt$$
And now that the definition for these to be uncorrelated, $$(f \star g) = 0$$  has to be true for all $\tau$. But what i don't understand is that unless either $f$ or $g$ are just zero for all $t$ there will be a value of $\tau$ where $$\overline{f(t)}g(t+\tau) \neq 0$$ It then follows that if there are any non-zero points then there must exist other ones which cancel so that the integral ends up being 0.

It seems to me that the chances of this happening for any useful finite signal is next to nothing and so it's highly unlikely that anybody could ever use the spectral subtraction.

Can anyone explain what I am missing?

Tl;DR version:  You are not missing anything; finite-duration signals cannot be uncorrelated signals.

If the crosscorrelation function $x\star y$ of $x$ and $y$ is zero everywhere, then the Fourier transform of $x\star y$, which is $X(f)Y^*(f)$ (or $X^*(f)Y(f)$ for left-handed folks), must also have value $0$ for all $f$. But, finite-duration signals have Fourier transforms with support $(-\infty, \infty)$
and so the product $X(f)Y^*(f)$ cannot be zero everywhere.

It is of course possible to have uncorrelated signals, but these have Fourier transforms with finite (and non-overlapping) support, and so (by duality) must have infinite duration.


Prepare IQ Data for Goertzel Algorithm I have a quadrature (IQ - 0..255) data byte stream that I receive from an rtl-sdr dongle as IQIQIQIQ.....

That means that negative values are < 128 and positive values >= 128

I am attempting to decode a 2FSK signal:

To normalize the data I do the following and get a resulting complex number with values ranging between -1 and 1 ... (more or less)

 Real = (I - 127.4) * 1/128
 Imag = (Q - 127.4) * 1/128

I demodulated with GNU Radio:

The scope plot correspond to the data that I get from a CC1101 tranceiver with an Arduino:
00000000011101110100010001110111010001000111011101000100010001110100010001000111011101000

This is the parameters for the CC1101 (SmartRF Studio)

The project is for a USB receiver to receive and identify signals from pager buttons. Unfortunately the CC1101 only works well for about 60% of the pagers.
That is why I want to use a RTL-SDR dongle and do my own demodulation.

GNURadio is great but it is a few hundered MB there is no minimal install ;-(

I would like to use the Goertzel algorithm to demodulate but it requires Real(float) input (not complex). If this is possible I could translate and compile the source to C or Delphi and create an executable that is maybe 1MB is size and a lot easier to deploy.

How do I prepare/use this data for the Goertzel algorithm?
Here is the prototype of the Goertzel C# code posted in another thread, that I want to use.

float goertzel(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)

It requires the data as an array of float. float* data
https://dsp.stackexchange.com/questions/23749/normalizing-the-magnitude-of-goertzel-filter-c

From reading the code I can see that it is not an array containing complex numbers (just an array of float) and I don't know how to get my IQ data into the format required by the Goertzel algorithm.

Finally my appologies and thanks to Marcus Müller (answered below) who spent a lot of time and effort to answer my original post in great detail, which had some incorrect, incomplete and confusing information, which I have now updated to be more complete.

Goertzel is defined on complex numbers $I + jQ$, not real numbers. So, directly apply Goertzel to your baseband.

Your 1/128 scaling is unnecessary, since that's only a factor and doesn't contribute to the decision.

Also, your -127.4 (which almost certainly is a misunderstanding of some sorts) just kills the bias of unsigned integers coming out the device. Since that all ends up in the (not evaluated) DC bin of your DFT, you don't have to do that, either, and don't gain anything by doing so, either (numerical accuracy of floating point operations aside, but with any reasonable observation length, i.e. any reasonable bit rate, your accuracy cannot be degraded - after all, you start with just 8 bit of precision).

So, take these numbers, and feed them into your Goertzel, just converted to the data type your Goertzel implementation needs.

Another note: don't let yourself be fooled, computers are fast. 1 MS/s is not very much, and you need at most a full symbol-length FFT, i.e. a $\frac{10^6}{2400}=416\frac23\approx 417$-point FFT (yes, there's FFTs for non-power of two) will also take nearly no time. I've just benchmarked this, on my really elderly AMD CPU:

I only got 60,000 transforms per second, corresponding to 60,000 baud – that's only 25 as fast as you would actually need, so that's OK, but a bit disappointing.

Now, note that FFTs work better, the smaller the prime factors of their length is. Under that point of view, 1 MS/s is a terrible choice – it's not a multiple of 2400 /s. So, instead, use a much smaller sampling rate that is a multiple of 2.4 kHz. What about 240 kHz?
That would give you 100-FFTs, and of these, I get 1.3 million per second (you'd only need 2400 per second to follow your bits!).

So, the question really is why you came to the conclusion you needed 1MS/s for your 2400 bd signal in the first place. Unless it's absolutely unreasonably wide-spaced in spectrum, that's a waste of sampling rate. Your Dongle supports a wide range of sampling rates, among these multiples of your nominal baud rate – use such to keep your effort in range, or resample.

So, generally, regarding the question how to preprocess the data:

Use something that gives you samples that you can work with. I do like GNU Radio (that's what I used to build above flow graph), and it comes with all the tools you need. Even a weak personal computer-style CPU (like the ARM in a Raspberry Pi) will have plenty enough horsepower so that you can do a full FFT. Don't reinvent the wheel, though, but use existing implementations.

Of course, for detection of one of two potential frequencies, a Goertzel is actually more appropriate then an FFT. But if you're clever, you realize that actually you can get away with two Goertzels, one of which covers every frequency above 0, one every frequency below 0. That's not really any more efficient than a 3-point DFT, of which you just ignore the DC bin. So, I'd argue:

take your samples; don't scale or shift them, because that's just cosmetics (DC offset doesn't matter to a "is there more energy left or right in the spectrum, neither does scaling).
Like anyone who strives for acceptable SNR, you'd filter to only consider the part of the spectrum where there might be your signal. Use a real-tapped (ie. spectrum-symmetrical) band pass filter to

suppress the DC component we didn't suppress, alongside with LO leakage
select your bands of interest from noise

Notice that above filter might optimally be (depends on whether you only consider an AWGN channel) the conjugate pulse shaping filter, shifted up with a cosine
Decimate the output of that filter. You can decimate very drastically here, just make sure if aliasing occurs, it always happens within the same band where the original energy came from. If you're not sure you can achieve that, only decimate to Nyquist rate.
Circularly shift in frequency (by multiplying with a $e^{j\frac{2\pi}{x}t} $) so that one of your tones now is at DC, and the other is not.
Your 2 Goertzels now collapse to the full 2-DFT, which really collapses to

First bin: Moving average of samples
Second bin: Difference of successive samples


Spreading symbols with sequence before transmission I  have been learning IEEE 802.11ad and the following processing algorithm

I was wondering what is the need for spreading with Golay sequence. I am guessing it has to with something like spread spectrum in CDMA applications.

Does anyone know why and whats the purpose of using such a block after modulation block?

Thanks

It is just the control channel that is spread.  It appears that it is spread primarily to make the signal more robust in the face of frequency selective fading, which is probably a major problem at 60 GHz.  It is easy to use spreading to get the robustness because the data rate of the control channel is low.  In other words, the spreading is fulfilling two roles: converting the low control channel data rate to the high OFDM symbol rate, and increasing the control channel's robustness.