Where do you recommend getting "standard" images for testing algorithms? I've noticed that there are a LOT of websites which claim to have a standardized repository of images (for image processing algorithms). However, I don't know which ones are actually trusted by the research community, and which are not. I know some file formats hint of untrustworthiness (e.g. never use a JPEG since its compression is lossy, and it is probably not the original image). Many of these repositories appear to have the same image, but their actual files do not match pixel-for-pixel.
So, what is your recommended repository for getting standard images for testing image-processing algorithms (e.g. images like cameraman, boat, airfield, peppers, pirate, etc.)

Research communities don't "trust"; they take a set of data and consider its properties (sample size, representativity, quality, systematic biases...) and then pick that. There's no "universally trusted" thing, as a dataset that's perfect for one purpose might be very wrong for another.
So, you'll have to work this from the other side up: pick a set of images that well represent what you want to do! Is it improving the quality of compressed photographs of cattle? Different dataset then you'd use to assess a denoising algorithm for architectural photography!
If you want to do that works on most photographic content, for example, you'd probably end up using the same dataset as the JPEG used to develop their compression standards (of the same name); not only has that been proven to actually be a useful set, but also it allows a fair comparison to established methods that use the same.


Is there two version of ulaw and alaw encoding in audio compression? I am new to encoding in audios and i have searched for ulaw and alaw and found there is conversion from audio to this encoding CCITT alaw/ulaw and as i have known alaw/ulaw belongs to ITUT group not CCITT so Is there two standards for conversion audio to  alaw/mulaw ?

and this conversion about CITT ulaw

Yes. You can choose the µ parameter freely in µ-law encoding. So, there's as many different µ-law encodings as there are possible µ; of course, only very few are standardized.
The same applies to the A in A-law.
You might want to read the wikipedia articles on µ-law and A-law, instead of just the G.711 article.


Digital low-pass filters; best practice application I have time domain strain gauge data from a structure which is subjected to low frequency, random loading. The signal / noise ratio is pretty good, but I would like to try to improve on it as much as possible.

From the FFT, I know that the dominant loading and response is below a frequency of around 1 Hz, but my signal has a constant level of noise up to the 10 Hz Nyquist frequency, above which any frequencies should have been removed by the low pass filter that was used with the measurement system - but I have very little knowledge of how this works. Although I can determine the dominant frequency response from my data in it's current form, I'm interested in processing the time domain data for fatigue calculations, and so any additional noise will introduce error.

I basically have two questions that I'd like to ask anyone with a good understanding of modern measurement and signal processing systems;

1) How is the low pass filter implemented in best practice for a measurement system like this? Is this done with the standard RC circuit prior to analogue to digital converter, or is this done after digital conversion with a processor and a suitable algorithm?

2) What would be the most suitable method to remove the noise present in my data above 1 Hz? I've had a bit of success with simple smoothing algorithms in the past (such as localised regression fitting; LOESS), but wondered whether there are more appropriate methods that I am not aware of, as smoothing seems a bit crude!

Any responses or pointers in the right direction would be really appreciated,

Many thanks,

Trev

I think a digital low-pass filter, applied after sampling, should be a good solution. One benefit of digital filtering is that you can obtain very sharp transition regions, and large rejection of out-of-band signals. The cost is relatively high computational complexity, in the sense that your processing code will execute many arithmetic operations.

Another benefit of digital filtering is that you can make use of a very large, well tested array of algorithms to design the filer, without necessarily understanding how they work. Obtaining the frequency response of a particular filter is very easy, so you can quickly evaluate if the filter meets your needs.

For illustration, this is how you would design and evaluate the filter in Matlab. I'll assume that the sampling frequency is 20 samples per second, and you want the cutoff frequency at 1 Hz.

% define frequencies:
fs = 20;
fn = fs/2;  % Nyquist frequency
fc = 1;     % cutoff frequency

% define the filter
b = fir1(100, fc/fn);  % filter order is 100; cutoff frequency is 1 Hz

% plot its frequency response
freqz(b);   % this command plots the filter's magnitude and phase responses.

The response obtained is:

You can see that the response is really nice: it is linear in the passband, and the rejection is around 60 dB to 80 dB out of band. Using this filter, though, requires around 100 multiplications and 100 additions per each processed sample.


Working with the DCT I am having a very hard time to implement the DCT algorithm. I have quite a few requirements like it has to work with NxN matrix or at least power of 2, it has to be 2D, it has to produce same output as FFTW fftwf_plan_r2r_2d(FFTW_REDFT10) it has to use real data and I need DCT I, II and III. And it has to be fast!

1) I could use FFT to compute DCT and wikipedia mentions this:

"One can also compute DCTs via FFTs combined with O(N) pre- and post-processing steps. In general, O(N log N) methods to compute DCTs are known as fast cosine transform (FCT) algorithms."

What are those 2 steps?

2) If I use Apples vDSP library is that even a good idea? since it needs an array twice as big as the DCT array (2N with radix2). And also vDSP is 1D only so I would end up doing it for each row and column to get a 2D result.

3) Is it possible to use an algorithm that works on 8x8 blocks and adapt it for my needs?

I am very new to learning about DSP so any help is appreciated! Thanks!

ps: Does anyone have any sample code for what I need? would help me a lot

While this is not much of an answer, here's a few comments while you wait for an expert to reply.

If this is one of your first DSP projects, it might be a bit ambitious just to dive straight into the fast implementations - you might want to start out with a naive implementation and then work your way up to the faster algorithms.

For the FFT-based approaches, you could check out how scipy does it:
http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.fftpack.dct.html

There's a lot of papers on the subject, but they might be a bit tricky to chew through. A couple of examples through some brief searches:

http://www.cin.ufpe.br/~vak/tg/papers/Chen-dct.pdf (1977 paper on a fast DCT)
http://arxiv.org/pdf/1203.3442.pdf (fast DCT-II)

It would probably be easier to understand from a book - I recall reading about the DCT and some fast algorithms with pretty elaborate explanations in "JPEG Still image data compression standard" by W. B. Pennebaker and J. L. Mitchell (I'm sure you can find loads of books on the DCT, specifically, though).

Also, this seems to be how they do it over at FFmpeg: http://www.ffmpeg.org/doxygen/2.1/dct_8c_source.html


help with bandwidth diagrams for digital transmission I am taking a undergrad network course that uses Data Communications and Networking 5th edition containing diagrams and legends which I cannot decipher.
Specifically the diagram on the right, has y-axis P for power, then z-axis f/N.  Does this mean the fraction f/N (frequency/bit rate) or does it mean f or N.  Also the start of chapter 4 explains that ratio r is data rate/signal rate or r = N/S.
In the legend it gives r=1 so re-arranging the equation, shouldn't S_ave = N/1?

Also, it explains that Manchester and Differential Manchester take up twice as much bandwidth because of the voltage change during the bit but comparing the right hand side diagram below with the respective right hand side diagram above, how do you determine that it is double?  Is it area under the curve?

The $x$-axis in the diagrams is normalized frequency (normalized by the data rate). Note that they define $S$, which is the maximum baud rate, and $S_{ave}$, which is the average baud rate. You get the maximum baud rate if your data is such that your signal is changing at a maximum rate. For NRZ-L or Manchester that would be a sequence of alternating $0$s and $1$s. However, on average the baud rate is smaller than that maximum baud rate. This is what is meant by $S_{ave}$, which is less than $S=N/r$.

For baseband signals such as these, the (effective) bandwidth is the maximum frequency above which there are no relevant frequency components. Since the required bandwidth is proportional to the baud rate, and since both Manchester and differential Manchester codes have $r=\frac12$ whereas NRZ-L and NRZ-I have $r=1$, the Manchester codes need approximately twice as much bandwidth as an NRZ scheme.


Why are the systems equivalent? Why are those two systems equivalent?

System 1:$$\boxed{x[n]}\longrightarrow\boxed{\uparrow 3}\longrightarrow\boxed{\downarrow 4} \longrightarrow\boxed{y[n]}$$
System 2:$$\boxed{x[n]}\longrightarrow\boxed{\downarrow 4}\longrightarrow\boxed{\uparrow 3} \longrightarrow\boxed{y[n]}$$

By drawing in time the first system's output my result for the first $4$ samples: $$y[0]=x[0], y[1]=0,y[2]=0,y[3]=x[12]$$

by drawing in time the second system's output my result for the first $4$ samples: $$y[0]=x[0], y[1]=0,y[2]=0,y[3]=x[4]$$

Why the statement is correct?

Upsampling (by $M$) and downsampling (by $N$) are interchangeable if and only if $M$ and $N$ are co-prime (i.e., they have no common factors other than $1$). This is the case in your example because $2$ and $3$ are co-prime.

Think about your specific example in the time domain. Introduce a signal $u[n]$ between the upsampler and the downsampler. In the first case you have

$$u[n]=\begin{cases}x[n/3],& n=3k,\quad k\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}$$

and

$$y[n]=u[4n]=\begin{cases}x[4n/3],& 4n=3k,\quad k\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}\tag{1}$$

In the second case you have

$$u[n]=x[4n]\tag{2}$$

and

$$y[n]=\begin{cases}u[n/3],& n=3k,\quad k\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}\tag{3}$$

With $(2)$ you get

$$y[n]=\begin{cases}x[4n/3],& 4n=3k,\quad k\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}\tag{4}$$

which is identical with $(1)$.


How to find the fundamental frequency of a discrete signal using partial autocorrelation? I hope you can help me with this question.

I am trying to calculate the fundamental frequency (to know the beats per minute) of a cardiac pulse signal using partial autocorrelation.
I use a 12-bit ADC and taking 4096 samples at Fs=512 Hz.

Heres what I'm trying:

First I do an offset of -2048 to align the signal to 0 (since the ADC gives me values from 0 to 4096).

Then I multiply it by the Hamming window wich I calculated as

$$ w[n] = 0.54 - 0.46 \cos \left(2 \pi \frac{n}{N} \right) $$

To get this

Then I do the partial autocorrelation like this.

$$ C_{xx}[n] = \frac{1}{N} \sum\limits_{m=0}^{N-1-n} x[m]x[m+n] $$

And finally I look for the maximum peak (that is not in the 0 position).

In this case I'm getting the maximum value at $n$ = 506.

But I don't know what to do with that index of the maximum of 506 if what I'm trying to get are the beats per minute.

Thanks you for your help.

Very simple: your sampling rate is 512Hz and your sample time is the inverse of this: about 1.953 ms. Multiply with 506 and you get the time length of the fundamental: 0.988s.

So you have one beat every 0.988s which is 1.0119Hz or 60.711 beats per minute.


Impulse response function from pole-zero graph I know that that the poles of this system are: $s_1 = -1 + 2j$ and $s_2 = -1 - 2j$. I think to solve this problem I would need to find the laplace of this plot, with the formula:
$$
W(s) = \frac{\rm zeros}{\rm poles}
$$
Then I would need to get the inverse Laplace. I am stuck on how to get the zero when it is not provided in the plot.

You don't need the zeros, you just need to understand the relation between poles and the corresponding time domain function. A real-valued pole $s_{\infty}=\sigma$ corresponds to an exponential function

$$\frac{k}{s-\sigma}\Longleftrightarrow ke^{\sigma t},\qquad t>0$$

A double real-valued pole corresponds to an exponential multiplied by $t$:

$$\frac{k}{(s-\sigma)^2}\Longleftrightarrow kte^{\sigma t},\qquad t>0$$

[You can look up how this generalizes to poles of higher order.]

A pair of complex conjugate poles $s_{\infty}=\sigma\pm j\omega$ corresponds to a sinusoid with frequency $\omega$ multiplied by an exponential:

$$ke^{\sigma t}\cos(\omega t+\phi),\qquad t>0$$

The scaling $k$ and the phase $\phi$ depend on the numerator of the Laplace transform, but here they are irrelevant.

So for the given pair of complex conjugate poles, we know that the corresponding time domain function must be an exponentially damped ($\sigma=-1<0$) sinusoid with frequency $\omega=2$.


Is an image a random variable or a random process? Is an image a random variable, where each pixel is a realization of the same random variable? Or is an image a collection of random variables (if an image is of size mxn, there are mxn distinct random variables)?

In the latter case, are the RVs (pixels) considered i.i.d. (even though they are not)?

While there are some obvious exception to it ("static" pattern on a television screen, the "dark frame" noise pattern of a camera), images are rarely generated by random processes. Declaring that an image is drawn from such or such distribution or generated by such or such random process is just a post-hoc modeling decision, and there is no "ground truth" to validate or invalidate this choice; other than the performance of the machine vision / image enhancement etc. method derived from the modeling decision.

So you can view a $m \times n$ image as a single random matrix (I assume this is what you mean by your first alternative - the image as a whole being considered as a single multidimensional random variable); or you can view it as a random field (a collection of random variables indexed by $|[1, M]| \times |[1, N]|$). I've encountered the random field view more often than the random matrix view.

When using a random-field approach, you can view each pixel as i.i.d ; or you can introduce dependencies between pixel values through the use of a Markov Random Field model. These are not the only option - you could very well consider a two-layer model where one first random process assigns a region index to each pixel of the image, and then where the value of each pixel is drawn from a distribution indexed by the region id! No approach is "better" than another. The more complex the model, the more "plausible" the images it will generate, but the more intractable the computations might become. It often helps, when using this kind of statistical approaches, to draw a few sample images from the chosen distribution/random process, and look at them to have a good grasp of what kind of assumptions you have built into your model.


Iris recognition using filters? I am implementing an iris recognition system in matlab and have performed the following steps according to an algorithm proposed in a paper

Initial image

Median filtering
For each row in the filtered image, assign a value of one to the pixels that have the minimum value in that row and the other pixels to zero
Again median filtering to remove white spots here and there.
States that: Get the row and column that has maximum
number of ones. Then find the center (xo, yo) and
the radii (rx, ry) as follows: xo = the index of this
column, yo = the index of this row, rx = the
number of ones in this row/2, ry = the number of
ones in this column/2

My question is how do i do step 5 in MATLAB?

Given a matrix A that consists of zeros and ones only (which you should have after step 4 of your algorithm, you could implement step 5 via the following:

% find column with the maximum number of ones. sum down the rows of the
% matrix and find the column with the largest sum
[ry,xo] = max(sum(A));
% find the corresponding row; do the same thing except sum across the
% matrix's columns
[rx,yo] = max(sum(A,2));
% the values "rx" and "ry" above contain the sums of the row and column
% (respectively) that have the maximum number of ones. divide them by
% 2 according to your algorithm in step 5
rx = rx / 2;
ry = ry / 2;

After the above (which is untested), then xo, yo, rx, and ry should be populated with the values that your algorithm indicated.


Additive but not homogeneous continuous system? I am trying to find an example of a continuous system that is homogeneous, but not additive.
So far, the only example I could find was an example from this page, which describes system as $y(t) = \overline{x(t)}$, where (as I understand it) $\overline{a}$ is the conjugate of a complex number $a$.

This, however, will work for a scaling constant $a$ that is complex, but not real, and I am wondering if it is possible to find a system $S$ which is additive, but not homogeneous, and this holds even when we scale the input signal with a real constant?

I have been thinking a lot about it and I could not come with an example, nor could I prove it does not exist

You are misinterpreting the issues entirely.

  A system with complex inputs and complex outputs that maps its input
  $x(t)$ to output $[x(t)]^*$, the complex
  conjugate of $x(t)$, is additive, is not homogeneous for complex scalars and is homogeneous for real scalars.

Consider the system with complex inputs and complex outputs that maps its
input to its complex conjugate.  Thus, it maps
$$x_1(t) \to [x_1(t)]^*, \quad x_2(t)\to [x_2(t)]^*.$$
So, what happens if the input is the sum $x_1(t)+x_2(t)$? well, the
rule is that the input gets mapped to its complex conjugate and so
$$x_1(t) + x_2(t) \to [x_1(t) + x_2(t)]^* = [x_1(t)]^* + [x_2(t)]^*.$$
Since this holds for arbitrary choice of $x_1(t)$ and $x_2(t)$, we
conclude that the system is indeed additive: the response to a
sum of signals is the sum of the responses.

Let $\alpha$ denote a complex
number. Then, this "complex conjugator" system maps
$$\alpha\cdot x(t) \to [\alpha\cdot x(t)]^* = [\alpha]^*\cdot [x(t)]^*
\neq \alpha\cdot [x(t)]^*$$
unless $\alpha$ happens to be a real number.  So, the system is not
homogeneous since the required condition does not hold for arbitrary
(complex) $\alpha$.
However, it is true that if $\alpha$ is restricted to be a real number,
then the scaling property (a.k.a. homogeneity) does hold.

Addendum:  The OP wants

  a system that maps strictly real inputs to strictly real outputs, and that cannot be scaled with a complex (only with real) constants, and that is additive, but not homogeneous. That's my initial question,....

Since the input and output are required to be real, and not complex, the
notion of "cannot be scaled with a complex constant" does not make sense
any more because the scaling of the input makes it a complex-valued signal,
and the system allows only for a real-valued input. So essentially the question
boils down to

"Does there exist a system with real-valued inputs and outputs
that has the additive property but not the homogeneity property?"

to which the answer is No, as pointed out in the
Wikipedia article on linear systems. Most of the ideas there are reproduced
in the dsp.SE answer
mentioned in a comment by MBaz but the dsp.SE answer goes on to make
the additional claim (not in Wikipedia) that

"not much else is needed to extend to $a \in \mathbb C$"

which is false.  Indeed, as MattL notes carefully in a comment,
MattL's answers are applicable to scaling by real numbers only.

An outline of the proof that for real-valued systems, additivity
implies homogeneity is as follows:

Consider a system with real-valued inputs and outputs that
maps any given $x$ to $A(x)$. The function $A(\cdot)$ has the
property that $A(x_1+x_2) = A(x_1)+A(x_2)$ and the system is called
additive since the response to the sum $x_1+x_2$ is the sum of
the responses $A(x_1)$ and $A(x_2)$.  Mathematical induction extends this to
$$A\left(\sum_{i=1}^n x_i \right) = \sum_{i=1}^n A(x_i).\tag{1}$$

Choosing all the $x_i$ in $(1)$ to have value $0$, we see that
$A(0) = nA(0)$ and so $A(0)=0$.
Choosing all the $x_i$ in $(1)$ to be $x$, we see that $A(nx)=nA(x)$
and so the scaling property holds for positive integers $n$.
Choosing $n=2$ and $x_2 = -x_1$ in $(1)$, we have that
$A(x_1+(-x_1)) = A(x_1)+A(-x_1)$. But $x_1+(-x_1)=0$ and so we
see that $A(0) = 0 = A(x_1)+A(-x_1)$, that is,
$$A(-x_1)=A((-1)x_1) = -A(x_1) =(-1)A(x_1).$$ So the scaling
property holds for $-1$, and by mathematical induction for
all negative integers as well.
Scaling holds for rational numbers $\displaystyle \frac 1n$ as well. Choose
$\displaystyle x_i = \frac 1n x$ in $(1)$ to get
$$A(x)=nA\left(\frac 1n x\right) \Rightarrow
A\left(\frac 1n x\right) = \frac 1n A(x).$$
Mathematical induction allows us to extend this to more general
rational numbers $\displaystyle \frac mn$ as well.
Finally, to show that scaling applies to arbitrary (not necessarily
rational) real numbers, Wikipedia says that we need to assume that
the mapping is continuous and use the fact the rational numbers
are a dense subset of the reals.

So, except for this last part which requires some mathematics
beyond high school, we can show straightforwardly that
additivity implies homogeneity.


What benefits can we get from the use of GFDM? Existing literatures say that the generalized frequency division multiplexing (GFDM) has lower out-of-band emission (OOBE) levels compared to orthogonal frequency division multiplexing (OFDM). However, this is actually untrue in a fair manner. In the previous papers, the OOBE levels of these two schemes have been compared under the assumption of the same number of subcarriers. But such an assumption is very unfair!

Obviously, the total number of subcarriers of OFDM system (denoted by 'K') is determined by considering the coherence time and bandwidth expected. Therefore, under the same (or similar) circumstances, it is fair to set the number of samples in the GFDM symbols to be equal to 'K', which means that MN=K (M: number of subsymbols in GFDM system, N: number of subcarriers in GFDM system). In short, comparing the OOBE levels under the condition of K=N is quite unfair because GFDM requires the longer coherence time than OFDM in that case, and two systems have different spectral efficiencies. Instead, assuming MN=K is fair.

However, when MN=K, there is no gain of OOBE of GFDM systems compared to the OFDM one because of circular filtering, and I have already checked it up through both mathematical analysis and computer simulations.

So, I want to know what benefits we can get from the use of GFDM. As I know, the improved PAPR performance can be expected. But I'm not sure it is meaningful because of high implementation complexity of GFDM transceiver. In my opinion, the OFDM with PAPR reduction techniques may be preferable.

How do you think about it? Am I missing something? If it is, can you please let me know what benefits we can obtain from GFDM?

Thank in advance

In a system where subcarrier crosstalk is not defined by how well the parties are synchronized, synchronization can might more easily be implemented. Whether that really compensates for the more complex filtering structure is up to research or even questionable.

Point is that you can actually, in math and simulation bring down OOBE with GFDM, not even considering different comparison bases for the number of subcarriers:

Consider this: The pulse shape of a single OFDM subcarrier is inherently a sinc. Hence, the first side lobes of the outermost carriers can only be $\approx13$ dB lower than the desired signal, and will continue periodically, whereas such a statement cannot be universally made for GFDM.


Analog Hilbert transformer I know the FIR approach, I have seen IIR, to, but I'd like to know if it's possible to implement a Hilbert transformer in analog domain, i.e. with integrators instead of delays. Is it possible? If yes, how?

[edit]
After some searching on my own, I stumbled upon an answer in sci.electronics.design - "Analog delay??" (I'm not posting the link as it seems it involves some character string and I don't know what influence will that have). A simple search will do.

Still, this involves a translation of a digital delay into its analog counterpart, clock, etc, no continuous-time values. Could this mean that there is no other way to do it (analog)?

Yes, it is possible to implement a Hilbert transform in the analog domain, although it is difficult to get the desired phase shift over a significant bandwidth when you take component uncertainties and imperfections into account. There is some discussion of the topic in this document.

I'm not sure why you think you would use an integrator instead of a delay. The analog counterpart of a discrete-time delay is a continuous-time delay. This can be implemented in various ways, although as noted before, achieving precise delays in analog electronics is quite difficult. Some methods that I've seen before for inducing delays include introducing capacitive delays (e.g. with an RC filter) or using surface acoustic wave (SAW) devices.


Easy Filter with FFT and Convolution I'am pretty new to this topic so i dont have much experience with DSP.
I want to filter(highpass) a WAVE file, i have already programmed the FFT, and now want to filter the FFT vector via convolution (which is a multiplication in frequency domain ?!).

1) Can anyone give me an easy equation for such a filter?

2) If i have an equation i compute it, also Fourier-Transform it, an multiply it with my FFT vector and inverse Transform it an i have applied the filter?!

Actually, yes. As Hilmar pointed out, it might be cheaper to use an ordinary filter, but that depends on the circumstances.

If you have two signals $A$ and $B$ with Fourier transform being $\mathcal{F}$, their convolution $A\ast B$ will be

\begin{equation}
A\ast B = \tilde{\mathcal{F}}(\mathcal{F}(A)\cdot\mathcal{F}(B))
\end{equation}

With $\tilde{\mathcal{F}}$ being the inverse transform. This is the well known "convolution lemma". There is not much of a theory to see. In the Fourier domain, you just replace convolution by (point-wise) multiplication. Your question already contained the answer, to me it just appears that you were not sure about the details.


Difference between Contrast and Intensity of an Image I wanted to generate Gray Scale wedge image of 10 Levels in MATLAB and then increase and decrease its Intensity. By high intensity, I mean mapping increasing the level of gray scale intensity in an image from lower value to higher values.

The examples I have tried were performed on a Gray Scale wedge image of 10 levels:

0   28   57   85  113  142  170  198  227  255

I was using the function imadjust. For increasing the intensity levels I used:

 imadjust(grayStepImage, [.35 .75], [0.7 0.8])

The output:

179  179  179  179  184  192  199  204  204  204

and the result in image form:

Similarly, for lowering the intensity I used

imadjust(grayStepImage, [.35 .75], [0.1 0.39])

Original and modified levels are:

0   28   57   85  113  142  170  198  227  255

26   26   26   26   43   64   84   99   99   99

and the result in image form:

Was I actually increasing/lowering intensity of the image? How do those operations relate to the case of increasing and decreasing the contrast of an image?

What I know about contrast is,

  Contrast is defined as the separation between the darkest and
  brightest areas of the image. Increase contrast and you increase the
  separation between dark and bright, making shadows darker and
  highlights brighter.

Does that mean, if I want to increase the contrast of an image I should increase the high intensity levels of an image to much higher and lower intensity level to much lower?

I am confused with contrast and the intensity as they seem opposite to each other.
Kindly help me out.

On intensity

Simply said, it's hard to talk about "intensity" of an image. Every pixel has its intensity (for greyscale images, they are usual allowed range is [0, 255]), but the concept of image intensity does not exist. If you are doing some kind of image analysis, you could be interested in a parameter describing image intensities, e.g. mean intensity (like @geometrikal said) or distribution of image levels (which is related to contrast).

On contrast

If you presume an idealized situation where a greyscale image contains only pixels of two intensity values, one for background pixels and one for object pixels, the contrast of the object would be the difference of those values. What it means on the displayed image --> the higher the contrast, the easier it is to spot (find, locate) the object on the image (for the human eye). Again, you would usually talk about contrast of an object in an image, it's hard to talk about contrast of an image where it is not possible to define an object.

As an example, let's look at several 1-D signals (you can produce images similar to your examples from these):

1)   0   0   0   0  20  20   0   0

2)   0   0   0   0 255 255   0   0

3) 200 200 200 200 255 255 200 200

4) 200 200 200 200   0   0 200 200

I'll assume that the two pixels different from the rest in each example are object pixels, and the rest are background. We have several cases here (and I could make more):

brighter object on dark background, all intensity values are low (i.e. the image is rather dark), contrast is low
bright object on dark background, most pixels have low intensity, contrast is high (it's a black-white image)
bright object on darker background, all intensity values are high (i.e. the image is bright), contrast is low (but higher then in the 1st example)
dark object on bright background, most pixels have high intensity, contrast is high (but not as high as in the 2nd example)

Here's another example supporting the fact that the contrast is usually an attribute describing an object, either in relation to itself or to the surrounding. Another important fact to notice from this example is that contrast generally describes the difference between the intensity levels, but the precise definition depends on the application (purpose for the measurement). Let me paraphrase the contrast definition used in hierarchical image segmentation from P. Soille, L. Najman: On morphological hierarchical representations for image processing and spatial data clustering as an example:

  The internal contrast of a connected component (object) corresponds to the largest intensity difference between two adjacent pixels belonging to this connected component. The external contrast is defined as the smallest intensity difference between a pixel of the considered connected component and an adjacent pixel not belonging to the considered component.

Of course, a different application could use a different measurement for contrast.

On imadjust

I don't usually work in MATLAB, but from the documentation, it is used to map intensity range (intensity levels in given range) of an input image, to an intensity range specified for the output image. The default, no-parameter call, will increase contrast of the image. But, you can get a brighter/darker image with imadjust with either increasing or decreasing the object contrast. Let me demonstrate this on my 3rd example:

imadjust(I3, [0.8; 1.0], [0.0; 1.0])

should output 0 0 0 0 255 255 0 0. You would get a white object on black background. In general, the image is darker (lower intensity), but contrast is higher (black-white).

imadjust(I3, [0.8; 1.0], [0.0; 0.2])

should output 0 0 0 0 51 51 0 0. Fairly dark object on a black background. The image is darker (lower intensity) and the contrast of the object did not change.

imadjust(I3, [0.8; 1.0], [0.85; 0.9])

should output 216 216 216 229 229 216 216. Very bright object on bright background. The image is brighter (higher intensity), and the contrast of the object decreased (object intensity levels are more similar to the intensity levels of the background).

You get the actual ranges the function is working with by multiplying them with the maximum gray level (255). For example, a range [0.2; 0.8] is actually intensity range [51; 204].

One more thing to take care of is that the function clips the values outside the first intensity range, and maps them to the new low if they're smaller or new high if they're larger then the range. All of my examples actually include this: the first range starts from 0.8 which maps to intensity 204, but the intensity of 200 from the input image is mapped to the output low in all the images.

So, it's actually just a simple scaling of image intensities (with cut-offs). Also, the default call to imadjust with only an image as an input parameter should increase contrast. I'd say that imadjust(I2) wouldn't do anything (there's a maximum contrast in my second example).

On contrast enhancement

Quote from P. Soille: Morphological Image Analysis:

  Image contrast enhancement refers to accentuation or sharpening of image features so as to make a graphic display more useful for visualization or analysis of the image by the human eye.

He also emphasizes: with enhanced contrast image analysis by visual (human) inspection is easier. In idealized images like in my examples, a computer wouldn't gain much in term of image analysis difficulty.

E.g. for object extraction by thresholding, it would just mean different thresholds should be used to extract the object from the background. But, a human examiner would spot the object much easier on high-contrast images. This all changes some in real images (objects don't have uniform gray levels), so contrast enhancement becomes useful in image analysis.

There are several methods of contrast enhancement:

point-based techniques, where the local neighborhood is not important. They are based on the analysis of grey levels through the whole image
neighborhood based techniques, where local neighborhood of the pixel is important. Examples are white and black top-hat operators, and toogle contrast operator.
transform based techniques, where the filtering is done on a transformed image before using the inverse transform (e.g. filtering in the frequency domain after Fourier transform).

On intensity adjustment

In term of greyscale image, the intensity of the pixel corresponds to it's brightness. The greater the intensity, the greater the brightness. This also means that increasing intensity can be viewed as brightening the image (while decreasing intensity can be viewed as darkening the image).

I would describe the process of uniformly brightening the image as increasing intensity while leaving the contrast unchanged in the whole image.
This actually means adding a constant value to all the pixels.

Now, as the pixel intensity values have a predefined range, typically [0..255], this means that the maximum intensity exists. This inevitably means that some pixels (that were different from each other before) will become white (i.e. intensity 255).

This "naturally" happens when you take a photo aimed towards something very bright -- sun or other light. Away from the light, you might see some details, but at the position of the light/sun and around it, you will get only white pixels meaning that the intensity (amount of light when the image was taken) "hit" its maximum (displayable/storable) value.

The only way you wouldn't lose details by this operation would be if the original pixel intensities belong to only a part of possible intensity range, e.g. if original pixels are in the range [10, 150], the image could be brightened by up to 105 intensity levels before you start to loose details.

As imgadjust is meant to preform intensity scaling, it can do much more than just brighten/darken image. If you wanted to emulate brightening effect with imadjust, you could write something like:

imadjust(I, [0.0; x], [1.0-x; 1.0])

where x is any number between 0 and 1 (e.g. x=0.5 would brighten a [0..255] image by 127).

That said, this is an overkill for such a simple operation. I'm sure matlab has a elementary operation that adds a scalar to all elements of a matrix, so you could just use that :)


Object detection versus object recognition I have a slight confusion differentiating between object recognition and object detection. Some people say object detection is a sub-topic of object recognition? Can someone clarify the the difference between these two topics?

To the best of my knowledge:

Object Recognition is responding to the question "What is the object in the image" Whereas, Object detection is answering the question "Where is that object"?

Hope someone can illustrate the difference by also generously providing an example for each.

You kind of answered your own question.

Object Recognition: which object is depicted in the image?

input: an image containing unknown object(s)

Possibly, the position of the object can be marked in the input, or the input might be only a clear image of (not-occluded) object.
output: position(s) and label(s) (names) of the objects in the image

The positions of objects are either acquired form the input, or determined based on the input image.

When labeling objects, there is usually a set of categories/labels which the system "knows" and between which the system can differentiate (e.g. object is either dog, car, horse, cow or bird).

Object detection: where is this object in the image?

input: a clear image of an object, or some kind of model of an object (e.g. duck) and an image (possibly) containing the object of interest
output: position, or a bounding box of the input object if it exists in the image (e.g. the duck is in the upper left corner of the image)


IIR coefficients that simply boost the signal Is there a way to get IIR coefficients that only amplify a signal, nothing more? I suppose some phase shifting would be ok as well as long as it's not too noticeable.

Background: I am using a MiniDSP board and I'm attempting to use the filter section as a level booster to drive the sound into a compressor downstream. The MiniDSP allows you to type in whatever filter coefficients you want in the EQ section.

I will assume that your IIR filter is an order-2 biquad. But the reasoning remains the same if the implementation is different.

Original filter : $H(z) = \frac{b_0+b_1\cdot z^{-1}+b\cdot z^{-2}}{1 + a_1\cdot z^{-1} + a_2\cdot z^{-2}}$

If you want to boost your signal by G

$H_G(z) = G\frac{b_0+b_1\cdot z^{-1}+b_2\cdot z^{-2}}{1 + a_1\cdot z^{-1} + a_2\cdot z^{-2}}$

simply distribute G to the numerator

$H_G(z) = \frac{G\cdot b_0+G\cdot b_1\cdot z^{-1}+G\cdot b_2\cdot z^{-2}}{1 + a_1\cdot z^{-1} + a_2\cdot z^{-2}}$

Simply multiply your b0-b2 coefficients by G.

In the case where you simply want to boost the signal without filering it at all. Set $b_0$ to $G$, $b_1$, $b_2$ to $0$, as well as  $a_1$ and $a_2$ to 0.


How to calculate the instantaneous frequency of a signal? I am trying to do the complex trace analysis of a seismic trace in Matlab where I need to calculate the instantaneous frequency of a discrete signal. I have used the formula according to Barnes. I am running an algorithm where I need to iteratively run complex trace analysis on the residual trace obtained after subtracting a matching wavelet from this trace. Though the formula gives me the results for lower frequencies <50 Hz, it does not give the right results for >60 Hz. I am not able to understand why this is happening. Can someone please tell me the reason behind this and how to do this correctly?

Instantaneous frequency is a poorly defined concept over wide-band (single channel) data, but can be a useful estimation given narrow-band signals of sufficient duration.  One trick for using an instantaneous frequency estimation algorithm which is usefully accurate only in one frequency band is to resample your data to move the frequency band of interest into the bandpass of the algorithm.  You may need to bandpass filter your data before resampling so that the source band ends up entirely within the algorithm's bandwidth.

In your case, you can try high-quality resampling your data to double the sample rate to extract instantaneous frequency estimates from 50 to 100 Hz.  Triple for 75 to 150, etc., assuming the sample rate and anti-aliasing mechanism captured spectrum in those ranges.


Image correlation confidence measure I am trying to calculate the disparity between stereo image pairs.  I'm using what appears to be the standard approach of correlating between a window of pixels from one image and a sliding window from the other.  The horizontal shift that gives the best correlation is the disparity.

However, the images contain large almost featureless areas.  In those areas, the calculated disparity is frequently wrong.  I would like to introduce a correlation confidence measure, so that the results in these areas can be discarded.

Is there a common measure of confidence that is suitable for this application?  At the moment I am struggling to distinguish between truly featureless areas, and those with low contrast that still have recognisable features.

First, on distinguishing low contrast features from noise, then on a possible alternate course of action to solve what I think you are trying to do.

From my experience, distinguishing between noise and low contrast features is a hard thing to do. The problem arises quite often when you try to do edge detection on an image. Using the most common algorithms such as Canny, Sobel, LoG, etc. generally results in a good deal of noise, which is dependant of the parameters you set (e.g.: a 'threshold' needs to be set (chosen) in order to define at which intensity value an edge-candidate pixel is classified as part of the edge or not, so what if the noise is more 'intense' than the correct edge?)

So, what's the point with stereo correspondence/disparity? Some scholars tackled on the problem of noisy edge detection and there is right now a pretty good edge detector, namely the Berkeley edge detector, which use machine learning techniques to 'improve' the edge detection.

Here are 3 images, left is original, centre is with Canny edge and right is with Berkeley. (Images credit: ETHZ-Shape-Class Dataset)

You can see that (while maybe some tuning could have been done on Canny) Berkeley is providing a 'better' (less noisy) edgemap than Canny.

So, you could see that there are techniques that you could possibly try to adapt in order to tell the salient features from the noise. I can't guaranty it will work since I haven't done it myself. Also, please note that these machine learning algorithms are pretty slow which prevent their use in realtime processing of images.

This said, what you really want to do is measure the dissimilarity between the to images by mean of correlation techniques, which was not successful apparently. Let me tell you about what is known as keypoint-descriptors. You could read on this topic there. The broad idea is to be able to evaluate the correspondence of some (a lot of) features between two images. It is often used in object detection (where you want to find the presence and location of a certain object in an input image) or in image stitching (such as panoramas maker, where you want to 'align' two or more images). Now, how can you tell that you have two images corresponding to the same object? These methods using keypoint-descriptors are, in the end, giving a measure of similarity between both images. Actually, lots of them are measuring the dissimilarity between the features to tell if they are similar or not...

I hope you will be able to better understand the situation with this. This is not a direct pre made solution, but it could help you in your work!


White noise and its phase after narrowband bandpass filtering Suppose $w(t)$ is a white noise with uniform spectral density $N_0/2$.  It is bandpass filtered to get a narrowband white noise $n(t)$ for $f_c - B_T/2 \leqslant |f| \leqslant f_c + B_T/2$. We write this in polar form by splitting into inphase and quad-phase wrt $2\pi f_ct$. $n(t) = n_I \cos(2\pi f_ct) - n_Q \sin(2\pi f_ct)$ such that we get a phase $\phi_n = \tan^{-1} (n_Q/n_I) $.
Is this phase uniformly distributed on $[ 0, 2\pi ]$? If yes, why?

...And now for a differing opinion....

The OP's representation of bandpass white noise as
$$n(t) = n_I \cos(2\pi f_ct) - n_Q \sin(2\pi f_ct)\tag{1}$$ is inadequate; because each sample path of this noise process is a pure sinusoid of fixed frequency $f_c$ Hz which is not noise-like at all. Why so? Well, a sample path is what one gets when all the random variables are assumed to have taken on appropriate values, and so with $n_I = \frac{\sqrt{3}}{2}$ and $n_Q = \frac 12$ (say), we get the sample path
$$\frac{\sqrt{3}}{2}\cos(2\pi f_ct) - \frac 12 \sin(2\pi f_ct) = \cos\left(2\pi f_ct+\frac{\pi}{6}\right)$$ which is a pure sinusoid as claimed.  A better representation of $n(t)$ is
$$n(t) = n_I(t) \cos(2\pi f_ct) - n_Q(t) \sin(2\pi f_ct)\tag{2}$$
where $\{n_I(t)\}$ and $\{n_Q(t)\}$ are low-pass uncorrelated white noise processes with identical power spectral densities that have constant value for $|f| \leq \frac{B_T}{2}$. Note that a typical sample path from $(2)$ is "sort of" sinusoidal-looking at frequency approximately $f_c$ Hz,  but the instantaneous amplitude and instantaneous phase (as well as instantaneous frequency (derivative of instantaneous phase) are "slowly varying" functions of time. Here, "slowly varying" means in comparison to the approximate time period $f_c^{-1}$ of each cycle of the sinusoid.

All this is fine and dandy but has nothing to do with the question being considered: the distribution of the instantaneous phase $\operatorname{atan2}(n_Q(t),n_I(t))$. The WSS theory tells us nothing about the distributions of the random variables $n_Q(t)$ and $n_I(t)$, except, of course, when $n(t)$ is a Gaussian process in which case $\{n_I(t)\}$ and $\{n_Q(t)\}$ are independent Gaussian processes and we know for sure that $\operatorname{atan2}(n_Q(t),n_I(t))$ is indeed uniformly distributed on $[0,2\pi)$. But, I contend that since we are discussing white noise processes, any distribution other than $\mathcal U(0,2\pi)$ for the phase would mean that the noise exhibits a preference for some phases and a dislike for some other phases, and this idea sticks in my craw as being totally unrepresentative of white noise. Indeed, even the answer by @StanleyPawlukiewicz that $\arctan(n_Q(t)/n_I(t))$ is uniformly distributed on $[-\pi/2,\pi/2]$ is disturbing because it suggests that it might be the case that $n_I(t)$ is always positive while $n_Q(t)$ has a more symmetric distribution. We really need $\operatorname{atan2}(n_Q(t),n_I(t))$ for defining the phase; $\arctan$ is inadequate for the job.

So, there you have it, folks, The phase is uniformly distributed on $[0,2\pi)$ because it offends my sensibilities to have any other distribution for the phase, and that's all there is to it; ymmv.


Power Relation in DSB-FC AM So, in my class I have to learn how to derive the power relation in a Double Side Band, Full Carrier Amplitude Modulated Wave However, I have hit a wall somewhere.

The basic equation is this :

$P_{total} = P_{carrier} + P_{lower-side-band} + P_{upper-side-band}$
Then,

$P_{carrier} = \frac{E_c^{2}}{R}$ where $E_c$ is amplitude of carrier.
Now , when converting to R.M.S, I divide the numerator by $\sqrt2$   and I get: $\frac{E_c^2}{2R}$

$P_{lower-side-band} = P_{upper-side-band} = \frac{(\frac{mE_c}{2})^2}{R} = \frac{m^2E_c^2}{4R}$

The problem is that the notes provided to me by my teacher state the numerator is $8R$ and I do not know why. Can someone please tell me if I have made a mistake somewhere.

Consider your carrier signal as $E_c\cos W_{ct}$ and your message signal is $E_m \cos W_{mt}$.

The double side band plus carrier signal is,

$ E_c\cos(W_{ct})+(E_c\cos(W_{ct})E_m\cos(W_{mt})) $

$ E_c\cos(W_{ct})+\frac{1}{2}(E_cE_m[\cos(W_{ct}+W_{mt})+\cos(W_{ct}-W_{mt})]) $

$ Ec*Cos(W_{ct})+\frac{1}{2}E_cE_m\cos(W_{ct}+W_{mt})+\frac{1}{2}*Ec*Em*Cos(W_{ct}-W_{mt}) $

Carrier being $E_c\cos(W_{ct})$

Upper Side band being $\frac{1}{2}E_cE_m\cos(W_{ct}+W_{mt})$ and

Lower Side band being $\frac{1}{2}E_cE_m\cos(W_{ct}-W_{mt})$.

For a signal $A\cos W_t$, the average power is $(\frac{A^2}{R}T)*(\int^T_0 \cos^2 W_t) = \frac{{A_c}^2}{2R}.$

Hence,

Average Power of Carrier is $\frac{{E_c}^2}{2R}$

Average Power of Upper Side band is $\frac{{E_cE_m}^2}{8R}$

Average Power of Upper Side band is $\frac{{Ec*Em}^2}{8R}$


Symbol synchronization with unknown symbol rate I'm learning about the early late gate symbol synchronization and understand that it needs to be given the symbol rate ahead of time.
How could a multi-symbol rate system use an early late gate?
It is my understanding that if the rate is known, the error is computed by going through the signal every samples-per-symbol and computing the error and eventually the timing offset converges.
But if the system has multiple symbol rates, how can a receiver know what the symbol rate without a control channel?
The one idea I have is to first use a symbol rate estimator, but the problem is that there will be some error from the estimator, then the error gets used in the symbol synchronization. Ideally, there is some synchronization method that don't need to know the symbol rate ahead of time or one that is resistant to small errors in the symbol rate estimate.
Edit
This diagram is what I'm looking at:

The sampler block is the one which needs to know the symbol rate since it is taking samples every $T_{sym}$ seconds. But if the symbol rate estimator is off, then even if the timing offset $\delta$ is estimated correctly, the sampler will take samples like $t=\delta$, $\delta+\hat{T}_{sym}$, $\delta+2\hat{T}_{sym}$, ... . Is there a way to detect that the sampler is off and refine the symbol rate estimate?

The one idea I have is to first use a symbol rate estimator, but the problem is that there will be some error from the estimator, then the error gets used in the symbol synchronization. Ideally, there is some synchronization method that don't need to know the symbol rate ahead of time or one that is resistant to small errors in the symbol rate estimate.

I mean, you need some tolerance anyways. No two oscillators in this universe are exactly the same frequency. And oscillator price pretty drastically scales with required long-term accuracy... so, you must assume that your receiver doesn't have perfect symbol rate knowledge, anyway.
Symbol timing recovery methods are typically resistant to small rate errors – look at the early/late gate, especially: You'll just see that you're a bit early more often than late if you're overestimating your symbol rate.
You just need to continuously run that.
This really isn't much different to phase recovery methods, which typically can also correct small frequency errors.

How could a multi-symbol rate system use an early late gate?

Uff, in such a system you'll often use a bit more complex synchronization methods, which take information from preambles, and, just as mentioned above, will keep tracking the symbol rate.


Comparing the loudness of two samples I have two wave files each containing a sample of a note. I wish to compare, using a program, each sample and present how loud the first sample is with respect to the second if each sample is played back to the listener with the same physical conditions (volume, distance from speaker etc.). For example, if the second sample is twice as loud as the first, I wish to present "2" to the user. I'm not interested in being very accurate. I simply want to give a reasonable estimation.

Lets say the first file contains only a 1kHz tone with an amplitude of 0.1 (0.1*sin(2*pi*1000*t)), while the second contains a 1kHz tone with an amplitude of 1.0 (1.0*sin(2*pi*1000*t)). So the difference is 10db correct? Does this mean that the second sample will be, roughly speaking, perceived as being twice as loud, with respect to the first?

I understand that phons are with respect to spl and that therefore if I have only the two wave files to work with, I can't talk in such terms since I've no idea what spl existed at the microphone when the recordings were made - or more importantly, the spl that exists at the listeners ears when they hear the recording. However, lets say I have a 3kHz tone with an amplitude of 1.0 (1.0*sin(2*pi*3000*t)) in my first wave file, while the second file contains a 10kHz tone with an amplitude of 1.0 (1.0*sin(2*pi*10000*t)). Isn't it fair to say that the first file will be perceived as being louder? Should it be fair, based on the graph to say that the 3kHz tone will be perceived as being roughly 4 times louder? Since in the graph you have to, roughly speaking, increase the 10kHz signal by 20db to get the same loudness as that of a 3kHz signal.

Again lets say that the first file contains a 3kHz tone with an amplitude of 0.1 (0.1*sin(2*pi*3000*t)), while the second contains a 10kHz tone with an amplitude of 1.0 (1.0*sin(2*pi*1000*t)). Is it fair to estimate that the difference in loudness is a factor of 0.5? That is, the first file will be perceived as being twice as loud as the second.

This is tricky. Loudness is a very complicated perceptual phenomenon and playback of sound in a room is a very complicated physical phenomenon. To do this right, you would need calibrated sound pressure signal measurements at the listener's ears and then run these through a suitable perceptual model (here is a good overview article: http://www.tcelectronic.com/media/1014016/skovenborg_2004_loudness_m.pdf)

For sine waves you could use the ISO 226 curves but with calibrated signal they are only correct for 1kHz. Let's say compare to sine waves at 200 Hz and one has ten times the amplitude: If playing back the first one 30 dB SPL the second one will be about 4 times louder. If you do this at 80 dB SPL the second will only be 3 times as loud at best. This gets even more complicated when multiple frequencies are involved.

Things get a lot more messy if you include the playback system as well. For pure sine waves the perceived loudness (and energy at the ears) can vary by as much as 30 dB (factor of 8 or so) just by moving a few centimeters. That's a basic effect of room acoustics. Furthermore most playback systems are fairly non flat (either by design or by sloppiness). Most cheap set of computer speakers will not provide a lot of output at 10 kHz.

Your best shot is probably to assume the signal are "normal" (music, speech, something broad band), apply an A-weighting, assume some average playback level and calculate a loudness difference from there. This is a gross simplification but depending on the application this may or may not be good enough.


Time and sample rate for QAM I am trying to understand how to get time as I already have sample rate. I want to draw a graph of sample rate vs time. I have information such as carrier frequency, channel bandwidth. Can that be used to get time.

So far, I have understood this: Samples per symbol and number of symbols for QAM

In order to get the time, I would need to do this:

T = 1/Samples per symbol

Would this be correct.

Consider you have a bandwidth of B. The sampling rate $F_s$ which is related to B using the relation $F_s > 2B$. Now the sample time is $\frac{1}{F_s}$. So if you have N samples  then total time is $\frac{N}{F_s}$

Is that what you are looking for?


Recursive DCT implementation I am trying to understand the recursive implementation of the DCT on the input signal (used for LMS filtering) according to the pictures below.

The pictures and formulas are taken from the paper: "Adaptive Inverse Control: A Signal Processing Approach", by Bernard Widrow and Eugene Walach

Block diagram of DCT structure:

Normal Matrix-Vector DCT-formula

Recursive DCT-formula:

What bothers me is how can this recursive DCT implementation yield the same result as the normal computation by using matrix-vector multiplication when the recursive formula uses the signal $x_{k-n}$?

But on the other hand the normal computation only uses signals up to $x_{k-n+1}$.


1st and 2nd Gaussian derivatives for edge detection I have a project about edge detection, I made research on internet about filters and read many articles but there is an article talked about using Sobel and Gaussian derivatives then it mentioned "edge" and "ridge" detection successively. I found the equation for 1st and 2nd Gaussian derivatives ($3\times 3$ and $\sigma= 0.8$) and applied the before mentioned method but I noticed blurred image with slightly enhanced edges.
I still don't get the idea nor the exact role of using different Gaussian derivatives together for better edge detection, and how can I enhance curvilinear structure detection to get cleaner edges. If someone can break it down it would be appreciated.

An edge remains a concept that is a bit complicated to define, as it may involve a certain level of interpretation. For a pixel-wise point of view, I consider that a potential edge breaks down into three main features: it is singular (non-continuous, non-differentiable) across one direction, and more regular (smooth) in the other direction, at a certain scale.
The singular part can be enhanced by the derivative part of the 2D filter; a first derivative should be extremal at the singularity, a second derivative could be zero. The regular part can be strengthened by the smoothing aspect the filter. Both should be combined at a certain scale, driven by the "size" of the filter (support and shape factor of the function, here $3\times 3$ and $\sigma= 0.8$) for you.
The proper choice of efficient sets of filters to detect edges and contours requires to carefully adapt to the nature of your image (resolution, etc.) and the type of noise it suffers.


How to choose a FFT algorithm? There are many algorithms that can calculate FFT. Even with Cooley–Tukey FFT algorithm, different radix can be used and the algorithms can divided into decimation in time and decimation in frequency. Any comment on how to choose these algorithms in practice?

Both Cooley-Tukey and Radix DIT & DIF are based on the same principle, dividing the N samples into two groups, and doing the same for the resulting two groups recursively. DIT and DIF generally use Radix2, that is, split N into two N/2 groups and provide a $N\log{N}$ time, while Cooley-Tukey is a generalization which splits it into $N_1$ and $N_2$ groups, where $N=N_1\cdot N_2$. Cooley-Tukey also provides $N\log{N}$ time in most cases.

Now, Radix DIT and DIF (both special cases of the Cooley-Tukey) are very similar and provide the same optimization in terms of operation savings, with the difference being bit reversed samples either on the input or the output. These can be done in $O(N)$ time, either before (DIT) or after (DIF) calculating the DFT.

Moreover, in some cases you don't even have do to the bit-reversal. Some applications (such as convolution) work equally well on bit-reversed data, so one can perform forward transforms, processing, and then inverse transforms all without bit reversal to produce final results in the natural order.

Finally, as an answer to your question, you can use Radix2 DIT or DIF FFT, depending on the application (see above paragraph), with the same result: Instead of:

$4N^2 \quad\text{real multiplications}\\
4N^2-2N \quad\text{real additions}.$

We achieved:

$2N\log_2{N}-7N+12 \quad\text{real multiplications}\\
3N\log_2{N}-3N+4 \quad  \text{real additions}.$

Please note that Radix DIT and DIF algorithms require that $N=2^p$, which is generally easy to satisfy. Increasing the radix gives us $\log_4$ for radix4, etc., which increases the calculation complexity on hardware, and in turn decreases the number of operations by ~$25\%$. Radix 2 and 4 are considered the most common, while Radix 8 (~$8\%$ optimization) and up generally demand too complex hardware for too small optimizations. Cooley-Tukey gives similar results, with more room for tweaking - in case of a specific $N$, for example.

For further optimizations, there are the split-radix algorithm which uses both Radix 2 and Radix 4 , and many more application-specific algorithms.

Update: Other algorithms include:

Bruun's algorithm:

It's not used much because Cooley-Tukey algorithm offers at least the same, and often better optimization. It was developed as a means to calculate DFT of real data (as opposed to complex).

Rader's algorithm:

Rader's algorithm is used for cases where $N$ is a prime number, often a big prime number. Note that above mentioned algorithms require $N = 2^k$, so they don't really work for other $N$. Rader's algorithm guarantees $N\log{N}$ time for prime $N$.

Bluestein's algorithm:

Bluesten's algorithm calculates the DFT via convolution, and is another way of dealing with $N \neq 2^k$. For prime numbers, it offers $N\log{N}$ time, for even numbers Cooley-Tukey algorithm is better, and for odd numbers it offers a somewhat better time than Cooley-Tukey.

Please note, however, that for all such cases for $N$ it's possible to zero-pad the input series until we get an $N = 2^k$ and then use the Cooley-Tukey's algorithm. Of course, for prime $N$ it's more efficient to use Rader's or Bluestein's algorithms.

Moreover, these algorithms were developed at different points in time, so some are more efficient than the others.


pyrtlsdr: Why is my output different than universal radio hacker? I'm trying to build a script to do some signal processing from my RTLSDR. But, the bits I'm getting is not what I would expect.

The preamble using 'Universal Radio Hacker' looks like this:

But after plotting all values my script is seeing (-1 until +1) it looks like this:

Often I'm also getting plots that don't make any sense. At least the last one has some resemblance of the URH preamble, but also a lot of plots don't.

What am I missing.. 'Universal radio hacker' also uses pyrtlsdr, so I know it can work.
- Complex, real, imaginary - Can't say that I understand it fully. Something wrong there ?
- Noise, gain.. Not 100% sure what to do there. URH seems to set gain to 25, but when I do that I don't get anything, so I left it at 'auto'
- I should be getting (almost) all of the samples. In the beginning I did lose some, but after I moved the postprocessing logic into separate processes it seems fine.

My script:

from rtlsdr import RtlSdr
from datetime import datetime
import asyncio
import matplotlib.pyplot as plt
from multiprocessing import Process, Queue

async def get_samples(sdr, q1):
    counter = 0
    timestamp1 = datetime.now()
    # Get this many samples every time
    async for samples in sdr.stream():
        q1.put(samples)
        counter += 1
        timestamp2 = datetime.now()
        if (timestamp2 - timestamp1).total_seconds() > 1:
            # To see if we missed any samples. Should be close to srate samples p/s
            print(datetime.now().strftime("%Y-%m-%d %H:%M:%S"), 'GETSAMPLES -', int((counter * 131072) / (timestamp2 - timestamp1).total_seconds()), 'samples p/s')
            counter = 0
            timestamp1 = datetime.now()

def check_samples(q1, q2, noise, samplesperbit):
    while True:
        values = []
        withinnoisecounter = 0
        signal = False
        if q1.qsize() > 10:
            print(datetime.now().strftime("%Y-%m-%d %H:%M:%S"), 'CHECKSAMPLES - Warning: Queue depth is', q1.qsize())
        samples = q1.get()
        for sample in samples:
            if not signal:
                # Start when signal breaks out of noise range
                if sample.real > noise or sample.real < -noise:
                    print(datetime.now().strftime("%Y-%m-%d %H:%M:%S"), 'CHECKSAMPLES - signal detected')
                    signal = True
                    values.append(sample.real)
            else:
                values.append(sample.real)
                if not (sample.real > noise or sample.real < -noise):
                    withinnoisecounter += 1
                elif withinnoisecounter > 0:
                    withinnoisecounter = 0
                if withinnoisecounter > samplesperbit * 3:
                    # Break when signal is within noise range 3 * samplesperbit
                    break
        if signal:
            q2.put(values)

def decode_fsk(q2, srate):
    # ToDo, now it just plots
    while True:
        values = q2.get()
        plt.figure(figsize=(20, 5))
        # Preamble should be in here
        plt.plot(values[:2000])
        plt.show()

def main():
    noise = 0.9
    # 38400 * 30 to get exactly 30 samples per bit.
    srate = 1152000
    samplesperbit = 1000000 / 38400 / (1000000 / srate)
    sdr = RtlSdr()
    # Just like in URH
    sdr.freq_correction = 1
    sdr.sample_rate = srate
    sdr.center_freq = 868.200e6
    sdr.gain = 'auto'
    # Run check_samples in another thread to make sure we don't miss any samples
    q1 = Queue()
    q2 = Queue()
    p1 = Process(target=check_samples, args=(q1, q2, noise, samplesperbit))
    p1.start()
    # Run decode_fsk in another thread to make sure we don't miss any samples
    p2 = Process(target=decode_fsk, args=(q2, srate))
    p2.start()
    # This is the main loop
    loop = asyncio.get_event_loop()
    loop.run_until_complete(get_samples(sdr, q1))

if __name__ == "__main__":
    main()

Thanks!

Your signal recording clearly shows that you have long streaks of 1.0 – that probably means you're clipping.

Your signal is thus broken. Make a new recording with less gain.


How do I find transfer functions from a state space representation? Suppose I have a MIMO system in state space representation, for example:

$A=\begin{bmatrix}
1 &2  &3 \\
 4&5  &6 \\
 7&8  &9
\end{bmatrix}$

$B=\begin{bmatrix}
2 &3 \\
 5& 7\\
9 & 1
\end{bmatrix}$

$C=\begin{bmatrix}
6 &7  &11
\end{bmatrix}$

$D=0$

I have used random number to fill these matrices. I am using Matlab. Now, suppose I want to find the transfer function from the input $u$ to the to an output $x_2$ for example, how is it possible to do this?

I know that I can from this matrix create the state space representation. So, suppose I want the state space representation of a plant, I would do this:

G = ss(A,B,C,D)

and if I want to get the transfer function from it I could do :

G = ss2tf(A,B,C,D)

and so from here I could plot the frequency response:

bode(G)

but now, suppose I want to obtain the transfer function from the disturbance to the output, or the transfer function from the input $u$ to the output, how can I do this ?

[EDIT] For example, how do I obtain the sensitivity function from a state space representation of a MIMO system? Or the control sensitivity function?

or the transfer function from the input  to the output, how can I do this ?

Is that not what your $G$ is?  If you want to find the transfer function from just one element of $u$ to the output, then either delete the columns of $B$ that don't pertain to that element of $u$ and get your transfer function, or just look at the column of the transfer function that matches the element you want.

  but now, suppose I want to obtain the transfer function from the disturbance to the output,

Then you would make a column for $B$ (or make a new $B$) that represents the effect that the disturbance has on the system, and extract the transfer function from that.

Edit

I neglected to include the actual math that Matlab is doing under the hood.  This is twice bad -- once because I did it, and twice because I really don't like people who just push the Matlab buttons without understanding what's actually going on.

If you have a system in state space representation:
$$\begin{split}x_k = A x_{k-1} + B u_k \\
y_k = C x_{k-1} + D u_k\end{split}$$ then you can take the $z$ transform:

$$\begin{split}X = A\ X(z) z^{-1} + B\ U(z) \\
Y = C\ X(z) z^{-1} + D\ U\end{split}$$
Then (leaving out steps, the grader will have words with me, but it's all linear algebra) you can solve for $X/U$:
$$\frac{X(z)}{U(z)} = C \left(I z - A\right)^{-1} B + D$$

Note that a nice thing about doing it this way is that if $B$ and $C$ happen to be single column and single row, then you get a nice, traditional scalar transfer function -- but if they're multi-dimensional, you get a very natural matrix representation for a transfer function that just drops naturally out of the math.


All pass systems- Poles and Zeros If an All Pass system has poles at $0.8\exp\left(-j \dfrac{\pi}{4}\right)$ and $0.8\exp\left(j \dfrac{\pi}{4}\right)$, There would be $2\pi$ jumps in the phase (principal value) at angular frequency $\omega= \dfrac{\pi}{4}$ and $\omega= \dfrac{7\pi}{4}$, but why would there be $2\pi$ jumps also at $\omega= \dfrac{3\pi}{4}$ and $\omega= \dfrac{5\pi}{4}$ (Figure 5.24 in Oppenheim Schafer and Buck "Discrete Time Signal Processing 2nd Edition)?

I understand that in all pass systems with real coefficient frequency response, a pole at $z= z_1$ would mean a pole at $z=z_1^*$ and zeros at $z=\dfrac{1}{z_1^*}$, $z= \dfrac{1}{z_1}$. With this reasoning a pole at $0.8\exp\left(-j \dfrac{\pi}{4}\right)$ would mean zeros at  $1.25\exp\left(j \dfrac{\pi}{4}\right)$ and $1.25\exp\left(-j \dfrac{\pi}{4}\right)$. But this still does not explain the $2\pi$ jumps in phase at $\dfrac{3\pi}{4}$ and $\dfrac{3\pi}{4}$.

Please Help. Thank you very much in advance.

There are no phase jumps. The discontinuities in the phase plot result from the computation of the principal value of the phase, which requires the phase to be in the range $[-\pi,\pi)$. But this is just a convention and has nothing to do with the physical system. Apart from issues with the principal value, the phase of allpass systems does not jump. Instead, the phase always decreases monotonically, because the group delay of an allpass is non-negative (Eq. (5.98) in Oppenheim's book). Also, if you look closely, the first "jump" is actually not at $\omega=\pi/4$ but a bit before, so it is not at all related to the pole angle.

'Real' phase jumps occur with systems whose frequency response has zeros. At these frequencies the phase jumps by a value of $\pi$. But obviously this can't happen with allpass systems because their frequency response is constant and doesn't have zeros. (Don't get confused: an allpass transfer function does have zeros in the complex $z$-plane, but not on the unit circle; this is why its frequency response, i.e. its transfer function evaluated on the unit circle, has no zeros).


Correct signal analysis flow (Wind Turbine vibration)? I have vibration data from WT (gearbox) sampled at 25.6 kHz. From a previous post I received a lot of good recommendation.

What I am trying to do:

Find early stage bearing faults by using an envelope analysis aka. an high frequency resonance technique.
This technique can:

  Magnify the low amplitude, high frequency signal by filtering out low frequency components. .. fault indication frequencies can be determined in the envelope spectrum through a spectrum analyzer

-Fault Detection Of Rolling Element Bearings by C. Liu

Current method:

Apply bandpass filter (low: 5000, High: 6000Hz)
Rectify resulting signal (absolute value AND squaring to remove aliasing)
Apply lowpass at 6000 Hz

-> Fault detection: Do this for many other 10 second samples over time to look for changes in kurtosis, FFT plot etc.

Zoomed in images of the signals processed:

Questions:

Should I upsample before or after rectifying?
Does this process look correct?
Should I do another LP filter to "clean" up the signal?
Should I do another LP filter at a lower cutoff frequency to "clean" up the signal?

Thanks in advance.

My approach resembles this one, from -Fault Detection Of Rolling Element Bearings by C. Liu

I think the question as is, is unanswerable, but here I am making an attempt of correct it with these comments as answer.

  Should I upsample before or after rectifying?

For the algorithm, if you are finally assessing against a threshold in magnitude, you should not care about down|upsampling. Do not add that.

  Does this process look correct?

Strictly speaking, you are taking the algorithm from your document. Since your document has 200 pages, this is someway hardly to say if you are taking correctly the approach taken by the author.

  Should I do another LP filter to "clean" up the signal?

This depend on the specific fault you want to capture for the bearing mechanics. From there, you can decide which frequencies to focus, and how to assess them.

  Should I do another LP filter at a lower cutoff frequency to "clean" up the signal?

Again, please check specifically how your algorithm is built. Once you have it implemented, we could study modifications or alternatives.


How to compute regions of matrix Let's say we have the following matrix:

0 0 0 0 0 0
0 1 0 0 0 0
1 1 0 0 0 0
1 1 0 0 0 0
0 0 0 0 0 0
0 0 0 1 1 1
0 0 0 0 1 1
0 0 0 0 0 0

I want to calculate the size of every connected region of 1s and also get the corresponding matrix elements of each region. The size should be computed as (width + heigth) / 2.

In the example the bottom left region would have size (2 + 3) / 2 and the bottom right region would have size (3 + 2) / 2.

Additional to the size of the regions I also want the corresponding entries of the regions, i.e. region 1 has entries (2,2), (3,1), (3,2) and so on.

How can this be done in Matlab?

I'm looking forward for the answers and wish everybody all the best.

Lets take the given matrix as
$$A=\left[\begin{array}{cccccc}&0 &0 &0 &0 &0 &0 &\\
&0 &\color{red}{1} &0 &0 &0 &0 \\
&\color{red}{1} &\color{red}{1} &0 &0 &0 &0 \\
&\color{red}{1} &\color{red}{1} &0 &0 &0 &0 \\
&0 &0 &0 &0 &0 &0 \\
&0 &0 &0 &\color{blue}{1} &\color{blue}{1} &\color{blue}{1} \\
&0 &0 &0 &0 &\color{blue}{1} &\color{blue}{1} \\
&0 &0 &0 &0 &0 &0\end{array}\right]$$

Then use the MATLAB function bwlabel to find the connected components in the matrix.

cc=bwlabel(A);

when you run this in MATLAB you will get cc as,

$$ cc = \left[\begin{array}{cccccc}
  &0&      0&      0&      0&      0&      0& \\
     &0&      \color{red}{1}&      0&      0&      0&      0& \\
     &\color{red}{1}&      \color{red}{1}&      0&      0&      0&      0& \\
     &\color{red}{1}&      \color{red}{1}&      0&      0&      0&      0& \\
     &0&      0&      0&      0&      0&      0& \\
     &0&      0&      0&      \color{blue}{2}&      \color{blue}{2}&      \color{blue}{2}& \\
     &0&      0&      0&      0&      \color{blue}{2}&      \color{blue}{2}& \\
     &0&      0&      0&      0&      0&      0& \end{array}\right]$$

1 indicate connected region 1, 2 indicates connected region 2, etc

If you want 2D pixel locations use, find function in cc

i.e.:

[r c]=find(cc==1) gives the row (r) and column (c) locations of the 1st connected region.

[r c]=find(cc==2) gives the row (r) and column (c) locations of the 2nd connected region, and so on.

Now you can get what ever you want..


Adaptive filter weight adjustment I have 3 sensor inputs: $a(t)$, $b(t)$ and $c(t)$. I want to design a filter such that the weighted linear combination of the three is always a constant. Kind of like:

$$w_1(t)a(t) + w_2(t)b(t) + w_3(t)c(t) = k$$

So from my undergrad modules I think I need a adaptive filters. I can perform training to find $k$ and filter weights at initialization.

May I know what is the best filter to implement in my case? $a(t)$, $b(t)$ and $c(t)$ are not independent.

Kelvin

What I was getting at in the comments above is that the linear system $w_1(t)a(t) + w_2(t)b(t) + w_3(t)c(t) = k$ has an infinite number of solutions, so you need to state some criterion that allows you to choose a unique solution. I think you have pointed out a constraint that is worth examining.

The idea: note that the linear equation I gave above is the equation for a plane in $\mathbb{R}^3$:

$$w_1(t)a(t) + w_2(t)b(t) + w_3(t)c(t) = k$$

This defines a plane in three-dimensional Euclidean space.
$\left[a(t), b(t), c(t) \right]$ is a vector that is normal to the plane.
The collection of all vectors in $\mathbb{R}^3$ of the form $\left[w_1(t),\ w_2(t),\ w_3(t)\right]$that satisfy the above equation are points on the plane.
Since you said that $w_1(t), w_2(t), w_3(t)$ don't change much over a short time duration, then you can assume that the point $\left[w_1(t+\Delta t),\ w_2(t+\Delta t),\ w_3(t+\Delta t)\right]$ should be geometrically close to the point at the previous time instant, $\left[w_1(t),\ w_2(t),\ w_3(t)\right]$.

So, the algorithm would look something like this:

Initialize your algorithm by finding $k$, which you said you can do.
Solve for the initial weights $w_1(t), w_2(t), w_3(t)$, which you said you can do.
On subsequent time steps, measure $a(t)$, $b(t)$, and $c(t)$. This defines the plane in $\mathbb{R}^3$ that the filter weight vector can possibly lie upon.
Find the point on the plane that is closest to the filter weights from the previous iteration. Use this point as the new vector of filter weights.
Repeat.

I'm not sure if this will give the desired effect or not (as your inquiry is light on details), but it has an intuitive geometric explanation. It might be worth a try.


Using Fast Fourier Transform to determine musical notes Hi guys i am doing a course in Digital Filters and Spectral analysis. We are given a coursework/homework, and I have absolutely no idea what to do with it. I come from Maths background, never done any signal processing before, and since I am new to the university I don't really have anyone to ask.

Problem:

In this exercise you are required to use spectral analysis techniques to determine the musical notes played within a short audio sample (with sampling frequency 44.1KHz). The sample will comprise a short sequence of 5 chords, each comprising 3 or 4 different musical notes played concurrently. Each note comprises a fundamental plus a series of harmonics at multiples of the fundamental frequency. All of the notes in this exercise belong to the 12 note equal temperament scale.

Here is some piece of code I scrambled so far:

[x,fs] = audioread('sample1_va18535.wav');
fs % fs is the sampling frequency usually 44.1 KHz
wavplay; % Play audio
N = 4410; % 0.1 seconds at 44.1KHz
N1 = 2205; % 0.05 seconds
n = N1+1:N1+N; % 0.05 – 0.15 seconds
xn = x(n,1); % Select left channel of short clip
wavplay; % Play short clip
t = n/fs; % Time index
plot(t,xn);
xlabel('Time (s)');
ylabel('Amplitude');
window = hamming(N); % Create window
wxn = xn .* window; % Apply window
Xk = fft(wxn); % DFT
k = 0:N-1;
f=k*fs/N; % Frequency in Hz
pause;
plot(f,abs(Xk));
ylabel('Magnitude');
xlabel('Frequency (Hz)');

Now I am pretty sure this is far from over ( I am not even sure if this code is correct).Would anyone be so kind and explain to me how can I generate those frequencies? I know I am required to create a vector of fundamental frequencies and ignore any harmonics.

Here is the file: https://ufile.io/lyge5

Unless your audio data consists of pure sinusoidal note sources, a raw FFT magnitude plot will give you all the harmonics or overtones of pitched musical sounds, as well as any fundamental pitch frequency spectrum (if any).  So, with most common polyphonic music recordings, pitch detection/estimation, if possible, usually involves a significant amount of post-processing of any FFT results.  Not just picking the largest magnitude peaks.

Either Cepstral processing (or computing the Cepstrum), or the Harmonic Product Spectrum algorithm, look for periodic sequences of overtones or harmonic peaks embedded in FFT spectra, which might indicate a source pitch candidate.  (e.g. Look for embedded trains of some number of magnitude peaks at evenly spaced frequency multiples.)  But, especially with polyphony and many common chord types, there will likely be some number of false positives that you will have to figure out how to discount or ignore.

Or, instead of using an FFT, you could also try looking at auto-correlation peaks (or other lag-based partial similarity measures), and see if these peaks correspond to the appropriate pitch range and degree of polyphony.


Matlab : How to plot the constellation diagram for this method? I am trying to generate symbols +1,-1 using a function by partitioning the interval of the function into two levels. The function is the Tent Map. This in a way achieves bpsk modulation.  I want to see the constellation diagram for bpsk . But I am getting error

Undefined variable "comm" or class "comm.RaisedCosineTransmitFilter".

Error in GenerateData (line 32)
txfilter = comm.RaisedCosineTransmitFilter('Shape','Normal', ...

So, I typed in >>

help comm.RaisedCosineTransmitFilter

and I believe that my Matlab does not have this class because this is the output of the help

comm.RaisedCosineTransmitFilter not found.

How can I view the constellation diagram?

This is the way I generate +1/-1 symbols for bpsk.

%generate bpsk
 N = 10000; %number of data points
x(1) = rand(1,1);
for i =2:N
x(i) = 2 * min(x(i-1), 1 - x(i-1)); %Tent map implementation
end
   modData = 2*(x>0.5)-1;  %generating symbols
  % modData holds the modulated symbols
  qpsk = comm.QPSKModulator;
  Rup = 16;
txfilter = comm.RaisedCosineTransmitFilter('Shape','Normal', ...
    'RolloffFactor',0.5, ...
    'FilterSpanInSymbols',10, ...
    'OutputSamplesPerSymbol',Rup);
  constDiagram = comm.ConstellationDiagram('SamplesPerSymbol',Rup, ...
    'SymbolsToDisplaySource','Property','SymbolsToDisplay',100);

txSig = txfilter( modData);
constDiagram(txSig)

Could somebody please show the correct way to obtain the constellation diagram. I am using Matlab2011b version on 64 bit i7 processor. Thank you.

The function you are using is from the Communication Systems Design Toolbox. You might have it not installed/licensed on your PC.

In case you just want to plot a constellation diagram, you can go with this code (Note, that it works for QAM constellations, not BPSK (since BPSK has only one domain)).

values = -7:2:7;  % Values in real and complex domain
[X, Y] = meshgrid(values, values);
QAM = X(:) + 1j*Y(:);
plot(QAM, 'x')

Actually, the QAM constellation contains all combinations of possible real and complex values. Hence, you need the cartesian product between real and complex parts. Taking this solution, you can plot the constellation diagram.


Strategy / Method for Implementation of the Fastest 1D Linear Convolution / Correlation I am after the method and its implementation of the fastest convolution of 1d signals:

$$ y \left[ n \right] = \left( x \ast h \right) \left[ n \right] $$

Where $ \ast $ is the linear convolution operator.

Let's assume that we have infinite storage or memory but we have finite time. What strategy or algorithm can we use for a faster discrete time convolution operation?

Could someone compare the known methods for implementing convolution: Direct Convolution, Frequency Domain Convolution / FFT Convolution and Overlap and Save / Overlap and Add methods? Is there something even faster?

I compared 3 implementations for Linear Convolution of 1D signals:

Direct - Using MATLAB's conv() funciton.
Overlap and Save - Implemented in MATLAB with tuned loop to prevent allocation and optimal choice of the DFT window.
Frequency Domain - Using MATLAB' fft() and proper padding to implement Linear Convolution using Circular Convolution.

For various lengths of the input signal and the kernel I got this result:

I only compared cases the signal is not shorter than the kernel hence the upper triangle (Dark Blue) is invalid.
For the lower triangle I coded the fastest method as following:

Bright Blue - Direct.
Green - Overlap and Save.
Yellow - Frequency Method.

So, there is nor practical reason to use Overlap and Save.
For kernels with length up to ~400 samples it is better to use Direct Method.
For longer kernels it is better to use Frequency Domain.

Now, the actual number of samples as the border between Frequency Domain and Direct will be different for different CPU's and Memory Bandwidth (One can run the script on his computer to see).
Yet as a guideline I'd say:

If you can hold all the data in Memory - Don't use Overlap and Save.
Unless your kernel is longer few hundreds of samples use direct.
For kernels with more than few hundreds samples use Frequency Domain.

The full MATLAB code is available on my StackExchange Signal Processing Q52760 GitHub Repository.


Resampling signal in MATLAB I read a wav file in MATLAB and it has a sampling rate of 16kHz. I want to pass it through a filter I designed which has a sampling rate of 44.1 kHz. So, I thought I would use resx=resample(x,44100,16000) first in order to make them have the same sampling rate. I can see that I will have more samples in time now, but what is the expected change of the fft of my signal? I thought that the fft of the resampled signal would be the fft of the initial signal (0-8kHz) stretched to (0-22.05kHz) but according to the MATLAB plots it is not. Is it logical?

To add to Famous Blue Raincoat's comment: the process of upsampling a signal aims to increase its sampling rate without affecting the signal otherwise. In your example, outputting x or resx through your sound card should produce exactly the same sound. Upsampling is useful when you need to match sampling frequencies between different systems without affecting the signals otherwise.

Likewise, downsampling aims to modify the signal as little as possible. All frequencies below the new Nyquist frequency will be left untouched.


Kiss fft does not work after giving it more than 32 samples I am trying to take data from an accelerometer and apply Kiss FFT to the samples. I'm using a Freescale Kinetis FRDM-K22F board. I want to use 64 samples, but when I run the program I get an error saying "kiss fft usage error: improper alloc" I started turning down the sample size and saw that the FFT does work with 32 samples, but giving it 33 samples the program just stops and returns no errors. Giving it any more samples gives similar results.

I played around with how I set up the FFT and followed a few websites and forum posts:

KissFFT output of kiss_fftr

http://digiphd.com/programming-reconstruction-fast-fourier-transform-real-signal-kiss-fft-libraries/

Kiss FFT on a dsPIC33

From what I can see, I haven't done anything different from what the above websites and forums have done. I've included my code below. Any help or advice is greatly appreciated.

void Sample_RUN()
{
    int size = 33;
    kiss_fft_scalar zero;
    memset(&zero,0,sizeof(zero));
    kiss_fft_cpx fft_in[size];
    kiss_fft_cpx fft_out[size];

    kiss_fftr_cfg fft = kiss_fftr_alloc(size*2 ,0 ,NULL,NULL);

    signed short samples[size];

    for (int i = 0; i < size; i++) {
        fft_in[i].r = zero;
        fft_in[i].i = zero;
        fft_out[i].r = zero;
        fft_out[i].i = zero;
    }

    printf("Data Collection Begins \r\n");

for(int j = 0; j < size; j++)
{
    for(;;)
    {
        dr_status = My_I2C_ReadByte(STATUS_REG);
        dr_status = (dr_status & 0x04);
        if (dr_status == 0x04)
        {
            //READING FROM ACCEL OUTPUT DATA REGISTERS
            AccelData[0] = My_I2C_ReadByte(OUT_X_MSB_REG);
            AccelData[1] = My_I2C_ReadByte(OUT_X_LSB_REG);
            AccelData[2] = My_I2C_ReadByte(OUT_Y_MSB_REG);
            AccelData[3] = My_I2C_ReadByte(OUT_Y_LSB_REG);
            AccelData[4] = My_I2C_ReadByte(OUT_Z_MSB_REG);
            AccelData[5] = My_I2C_ReadByte(OUT_Z_LSB_REG);

            // 14-bit accelerometer data
            Xout_Accel_14_bit = ((signed short) (AccelData[0]<<8 | AccelData[1])) >> 2;     // Compute 16-bit X-axis acceleration output value
            Yout_Accel_14_bit = ((signed short) (AccelData[2]<<8 | AccelData[3])) >> 2;     // Compute 16-bit Y-axis acceleration output value
            Zout_Accel_14_bit = ((signed short) (AccelData[4]<<8 | AccelData[5])) >> 2;     // Compute 16-bit Z-axis acceleration output value

            mag_accel = sqrt(pow(Xout_Accel_14_bit, 2) + pow(Yout_Accel_14_bit, 2) + pow(Zout_Accel_14_bit, 2) );

            printf("%d \r\n", mag_accel);

            samples[j] = mag_accel;

            break;

        } // end if
    } // end infinite for
} // end for

for (int j = 0; j < size; j++)
{
     fft_in[j].r = samples[j];
     fft_in[j].i = zero;
     fft_out[j].r = zero;
     fft_out[j].i = zero;
}

printf("Executing FFT\r\n");
kiss_fftr(fft, (kiss_fft_scalar*) fft_in,  fft_out);

printf("Printing FFT Outputs\r\n");
for(int j = 0; j < size; j++)
{
    printf("%d \r\n", fft_out[j].r);
}

kiss_fft_cleanup();
free(fft);
} // end Sample_RUN

It is possible that the program expects a power of 2 samples. Does it work with 64 samples? The FFT by definition usually refers to power of two series.


Normalizing a vector with Mel Cepstrum Coefficients, Delta Coefficients and Log Energies for both This is my very first signal processing exercise so I'm very new to this space. That being said the goal is to do isolate voice recognition in identifying someone saying yes and no.As the title suggests I have a 1x26 vector that looks as follows:

  Columns 1 through 10

  (P)160.2504   66.7557  -35.9298  144.1598  -11.9133  -16.2184  -26.4730   78.0673   10.2075  -39.4318

  Columns 11 through 20

  -46.6702    5.9283  -42.9540    (P)0.2614    0.0382   -1.3315    0.1326   -1.4256    0.1936    0.4081

  Columns 21 through 26

   -1.4431   -1.6916   -1.1668   -0.2901   -0.4499   -0.4309

Where (P) denotes the log energy for the Mel Cepstrum Coefficients and Delta Coefficients respectively. This being said I would like to know how to normalize these values but how? I'm not sure what exactly to divide the energies by, let alone the coefficients. Like I said first project so a really go post guiding me would be REALLY appreciated!

You could collect a large collection of voice signals, extract the MFCCs, and compute the mean and standard deviation of each coefficient. Keep these numbers and use them once and for all.

Then, you can normalize the sequence of vectors from an utterance by subtracting the mean and dividing by the standard deviation you have pre-computed (do not recompute the mean and standard deviation on each utterance - since each file utterance will have slightly different statistics).

This might be useful if you use the euclidean distance to compare MFCCs - for example when using DTW. However, more advanced statistical modelling techniques do not need this normalization step - a gaussian model, for example, will actually learn the mean and covariance matrix of the coefficients without having to pre-normalize the data (Euclidian distance on mean/std normalized data is just like Mahalanobis distance with a diagonal covariance matrix).

For sake of accuracy, I should point out that this normalization scheme has a flaw... The last MFCC coefficients are known to be very noisy and less informative - this follows from the energy compaction property of the DCT (the less relevant stuff is in the highest coefficient), or from the view of the DCT as an approximation of a PCA on the log mel spectra. Thus, the mean/standard-deviation scheme I have presented might "blow up" the importance of the last coefficients. Machine learning/feature selection algorithms might be robust to this and would discard them anyway, but the euclidian distance is not. A compromise I have used in past research is to normalize by a power < 1 of the standard deviation ; or to scale down, after normalization, the nth MFCC coefficient by $\alpha^n$ where $\alpha = 0.95$. Other normalization schemes implement "liftering" - by multiplying the sequence of MFCC point-wise by a window function.


Coefficients of a Zero Forcing linear equalizer for infinite length input Given a discrete-time white noise model:
$$y_k = \sum_{i=0}^{\infty} \alpha^{-i}a_{k-i} + \eta_k$$
where $E(|a_k|^2)=2$, $|\alpha|>1$, $\eta_k$ are uncorrelated zero mean random variables with $E(|\eta_k|^2) = 2N_0$.

Find the Z-domain transfer function for this channel model.

Find the coefficients of a Zero Forcing linear equalizer for this system, that operates on the samples $y_k$ and removes all the intersymbol interference.

For what values of $\alpha$ the channel becomes free of intersymbol interference?

My attempt:

For part 1, $F(z) = z + \alpha^{-1}z^{-1} + \alpha^{-2}z^{-2} + \ldots$. However I don't know how to do the rest.

The following are hints as this is reads as a homework question. Please note that these are not ground breaking hints and can be found in a good text or summary on zero forcing somewhere I'm sure.
Hint 1
The first step would be to look at the model and be able to pick out what the channel is (you have channel convolved with signal plus noise). It seems you have identified the channel as, $h_k = \alpha^{-k}u_k$, where $u_k$ is the unit step function, and the z-transform can be found by directly applying the formula.
\begin{align}
H(z) &= \sum_{k=-\infty}^{\infty}h_k z^{-k} \\
&= \sum_{k=-\infty}^{\infty} \big(\alpha^{-k}u_k \big)z^{-k} \\
&= \sum_{k=0}^{\infty} \alpha^{-k} z^{-k} \\
&= \alpha^0z^0+\alpha^{-1}z^{-1}+\alpha^{-2}z^{-2}+\cdots
\end{align}
Just be careful when simplifying...
Hint 2
Zero forcing (ZF) means that the equalizer is trying to completely eliminate intersymbol interference (ISI), and ZF does achieve this in the noise free case. For the noise free case, in the z-domain your model becomes:
$$ Y(z) = H(z)A(z) $$
The ZF equalizer should be $W_{\text{ZF}}(z)=\frac{1}{H(z)}$ so that the equalized signal will have eliminated the ISI caused by the channel.
Hint 3
This should become clear when you consider what happens to each symbol as it goes through the channel. Take a single symbol being sent over the channel at $k=0$. The received noise free version is $a_0+\alpha^{-1}a_{-1}+\cdots$, or in words this is, $\text{Original symbol } + \text{ Previous symbol } + \cdots$. Can you pick an $\alpha$ to ensure the only non-zero value is the original symbol?


Interpolate/Decimate with a single filter? This is not a homework question (I'm out of school now 2 years). I'm thinking, let's say you have 2 systems, one at 32khz and the other at 48khz and you want to go between the two. Is there a way to accomplish this with a single filter? I have a feeling that if you make an FIR where every other coefficient is zero and exploit the fact that if fs = 32khz and 48-32=16=fs/2, then zeroing out every other coefficient of an FIR (but how many taps) would do both the zero-padding and the anti-aliasing simultaneously. Am I on the right track or is there some other trickery involved?

EDIT:

Now that I think about it, to go from 32khz to 48khz you interpolate 3x and decimate 2x (32*3=96khz...96/2=48). So in effect you would pad two zeros first but then take out one in the decimation, but just inserting one zero would yield 64khz decimation (with a 32khz lowpass decimation).

EDIT 2:

Furthermore, wikipedia interpolation:
"Upsampling requires a lowpass filter after increasing the data rate, and downsampling requires a lowpass filter before decimation. Therefore, both operations can be accomplished by a single filter with the lower of the two cutoff frequencies. For the L > M case, the interpolation filter cutoff,  \tfrac{0.5}{L} cycles per intermediate sample, is the lower frequency."

Yes, you can interpolate and decimate at the same time.  This is called "resampling".  If you google resampling you will find lots of information about it.  And yes, your reasoning about resampling is mostly correct.

When thinking about resampling theoretically you usually put the interpolation first to avoid Nyquist issues.  An interpolator is upsampling followed by a low-pass filter to get rid of the aliases.  A decimator is a low-pass filter followed by dropping samples.  You can combine the two by combining the two low-pass filters into one filter.  If the overall rate change is greater than one, the low-pass filter has the same pass-band/cutoff as the interpolation filter.  If the overall rate change is less than one then it has the same pass-band/cutoff as the decimation filter.  The filter operates at the interpolation sample rate.

You can make the filter more efficient by not actually upsampling the data or dropping filter outputs.  Instead, you use a polyphase filter for the interpolation and simply not calculate the results that will be dropped, but that is a topic for another day and another question.


Downsampling by non-integer factor Filter requirement I am trying to implement a logic to downsample data in non-integer factor.Lets say,converting 48Ksps to 1344 samples per second.Have two step implementation of upsample by K ,and downsample by L is in my mind, such that K/L gives me the required ratio.

The thing that troubles me is the filter(if any) required in between the upsampling and downsampling. Can anyone help me understand the necessity of this filter?
Also is there any better methods to do so?

If the data is already bandlimited to below half the new sample rate, then a low pass filter may not be required.  If not, low pass filtering is required to prevent aliasing, which can corrupt your result.  Low pass filtering can be combined with Sinc kernel interpolation of all the new points needed at the new sample rate, e.g. only one step, no up/down pair of sample computations required.


What is the value of last sample of a discrete sinewave I have a simple question. If I compute N sample values of a single cycle of a discrete sinewave (whose initial phase is zero radians), will the value of the last sample be equal to zero?

No.  Assuming the period of the signal is N samples since your indexing is zero based the $Nth$ sample will be indexed as $N-1$.

Your signal would be modelled by:

$$ x[n] = A \sin \left( \frac{2\pi}{N}n \right) $$

Thus $x[N]=0$ and $x[N-1] = x[-1] = -x[1]$.