What is the basic idea behind Fourier transform? What is the basic idea behind (discrete and continuous) Fourier transform (FT)? In short, what is the difference between discrete and continuous FT?

I have read multiple answers on the web related to this topic, but they are all long, provide long examples or cryptic equations. I would prefer not to read another long answer or an answer with the example of the circles or ellipses. I do not care about that. Ideally, I would like to read short answers which give the idea behind this concept. Your target audience should be high-school students (or, definitely, not math or physics college students), or, in general, people that do not have a very solid mathematical background and do not want to read equations. Hence, I hope it is clear that equations should be avoided as much as possible, unless they can be understood easily by high-school students.

I think the discrete Fourier transform (DFT) is typically the easiest to use as a base of explanation:

A DFT is a linear transform that operates on a time domain signal. After the transform, the signal is in “the frequency domain”, i.e the signal has now been broken down into a vector/signal which contains information about what sinusoids are contained in the signal, and how strong they are. That’s pretty much it, but I’ll continue on to be more thorough.

So let’s say I have a time domain signal that has a 200Hz size wave and some noise: assuming my sampling rate is sufficient, after I do my DFT, I’ll easy be able to see the sine wave as a peak in the frequency domain, whereas the noise won’t really “add up” to much in the frequency domain.

Of course it’s important to mention that once you do a DFT, there is also an inverse DFT operation that “unapplies” the transform; doing this returns you to the time domain and gives you your original signal back.

For a continuous Fourier transform, the idea is essentially the same, except it is just expanded to calculus rather than discrete mathematics/linear algebra

EDIT:

Okay, we'll try to make this even simpler using some pictures via MATLAB. I've added a good amount of words to explain in detail, but there is almost zero math in the text; hopefully that will be sufficient for you.

We'll start off by showing two signals: one is a sinusoid (with a frequency of 50Hz, that is it oscillates 50 times per second. If you need more information than that, please consult a high school text book), the other is that same sinusoid, but with noise added.

Since you asked what noise was, you can think of noise as random numbers added into the signal which "corrupts" the signal, i.e. it is no longer a pure sinusoid.  A plot of such signals in the  time domain looks like this:

On the top plot, we can see that there is clearly a sinusoid present.  In the bottom plot, it pretty much just looks like a bunch of random numbers.  A DFT essentially breaks a signal down into sinusoidal components. It does this by "looking" for the each sinusoid in the signal (it's a bit more complicated than that, but since you wanted limited to no math we'll skip it).

So now that we know that a DFT can "look" for signals, we take the DFT of our two signals, and plot the results to find the following figure:

In the top plot, we have the DFT of our pure sinusoid.  Recalling that I assigned it a frequency of 50Hz, we would expect there to be a peak at 50Hz, and indeed there is.  You may notice there's a negative frequency component.  There is more mathematics required to explain that which involves unit circles, which you said in your original question that you "don't care about", so we'll just leave  that be.

Now, on the bottom plot, we have the DFT of the sinusoid plus our random noise.  We can see that those two peaks are still there at -50Hz and +50Hz.  When we just looked at the time domain signal, it looked like "nothing" was there, just a bunch of random numbers.  Through using the DFT, we were able to show that there is actually a sinusoid buried in that noise, and we were able to identify it easily in the frequency domain.  This is a very basic motivation for how/why one might want to use a DFT.

Naturally, the inverse DFT (IDFT) does this operation in reverse, i.e. it maps a signal in the frequency domain to the time domain. Knowing that, it's easy to see that if we take the DFT of a signal, perform no other operations, and then eventually take the IDFT of the output, we are given back the  signal we started with.  This is because the DFT is what is called a linear transform.  If you'd like to learn more about that, you'll have to learn more about linear algebra.

This is just one application of the DFT, there is a myriad of other applications, and it is a VERY important subject/tool in signal processing.  The Wikipedia article for the DFT has a list of a wide range of applications. If you're interested in learning more, I'd start there.

If someone else would like to provide details for the continuous time Fourier Transform, please do! I think you'll find its more or less going to be the same, but a bit more mathematical in its explanation.

For completeness, here's the MATLAB code for the figures I've shown:

f = 50;               % sinusoid frequency
fs = 500;             % sample rate
t = linspace(0,1,fs); % create a time vector
% create a sinusoid
x = sin(2*pi*f*t);
% create a noise signal
n = 1.5*randn(size(x));
% add in the noise to the signal x
xn = x+n;
% create a frequency vector that we'll use to label the plot
f = linspace(-fs/2, fs/2, numel(t));

% create a some plots
figure(1);
subplot(2,1,1);
plot(t,x);
xlabel('Time (s)'); ylabel('Magnitude');
title('Pure Sinusoid (No Noise)');

subplot(2,1,2);
plot(t,xn);
xlabel('Time (s)'); ylabel('Magnitude');
title('Sinusoid Buried in Noise');

figure(2);
subplot(2,1,1);
plot(f,20*log10(abs(fftshift(fft(x)))))
xlabel('Frequency'); ylabel('Magnitude (dB)');
xlim([-250,250]); ylim([0 60]);
title('DFT of Pure Sinusoid');

subplot(2,1,2);
plot(f,20*log10(abs(fftshift(fft(xn)))))
xlabel('Frequency'); ylabel('Magnitude (dB)');
xlim([-250,250]); ylim([0 60]);
title('DFT of Sinusoid Buried in Noise');


How to sketch a frequency spectrum for an AM signal? So I have an AM signal, given by:
$$v(t)=50\cos(\pi 10^{6} t )+ 20\sin(\pi 10^3t)\cos(\pi 10^6 t)$$
I was asked to sketch the spectrum of the signal obtained, so I found the fourier transform of it, and that's what I got,

$$V(f)=25 \left[ \delta(f - 500k)+ \delta(f+500k)\right]-5i\left[\delta(f-999.5k)-\delta(f+999.5k) \right]-5i\left[\delta(f-500.5k)-\delta(f+500.5k) \right]$$

And so far I only knew how to sketch the first term, as shown below,

But I don't know how to sketch the terms with the $-5i$.

Since a signal's spectrum is complex, it can't be drawn in a single two-dimensional plot. What is usually done is to do two plots; one is the magnitude spectrum and the other is the phase spectrum.

When looking at signals, most of the time, the magnitude spectrum contains all the information that is needed and the phase spectrum is ignored. When looking at filters, the phase spectrum is important because one usually wants the phase response to be linear in the filter's passband.

So, what you need is the magnitude of $V(f)$:

$$|V(f)| = 25 \left[ \delta(f - 500k)+ \delta(f+500k)\right] + 5\left[\delta(f-999.5k) + \delta(f+999.5k) \right] + 5\left[\delta(f-500.5k)+\delta(f+500.5k) \right].$$

The magnitude spectrum is real, so it can be drawn in a regular 2-D plot. Note that all values of the magnitude spectrum are positive.


DC and Nyquist values are higher than all other bins I am using a mean centered signal and applying fft to it, however the dc and nyquist values are a lot higher than all other values. Wondering if I am doing anything wrong or need to add anything else. The first image is the signal and the other is the fft output.

load Dog_1_interictal_segment_0001.mat

%interictal_segement is a struct with a 16 x 10000 matrix named data
data=interictal_segment_1.data;

%take first two seconds of the first eeg channel
data=data(1,1:10000);
freq=5000;
N=freq*2;
hz=linspace(0,freq/2,N);

averaged_mean=data-mean(data);
channel_frequency=abs(fft(averaged_mean)/N).^2;

figure(1)
plot(averaged_mean)

figure(2)
plot(hz,channel_frequency,'ko-','markerfacecolor','m','linew',1)
xlabel('Frequency (norm.)'), ylabel('Power');

For a mean centered (DC removed, pure AC) signal the first bin of the FFT output should be an exact zero. In your plot, this zero output sample is possibly obscured by the squeezing of the graph display; i.e., if you zoom into the left side of your plot, you should see that zero valued first bin sample...

Furthermore, for a typical FFT computation of real data, the second half of the FFT bins, is a conjugate symmetric mirror of the first half (excluding the DC and / or the Nyquist bins). Note that the Nyquist bin is (about) at the center of the FFT output, and to the right of this Nyquist bin is the mirrored content. Note also that exact Nyquist frequency bin exists only for even length FFTs and not for odd length FFTs.

Looking into your output, the last and the first bins seem like mirrors; but that's wrong as I have indicated, your fist FFT bin must be a zero, and even if it's not zero, it won't be mirrored to the last bin. What's mirrored to the last bin is actually the second bin of the FFT output.


2D Fourier Transform Is it possible to get a 2D Fourier Transform by first taking 1D Fourier transform in first dimension and then performing another 1D Fourier Transform in the other dimension? If yes, can you explain

https://www.ques10.com/p/5768/state-and-prove-the-property-of-kernel-separatin-1/
It is possible due to the separability property of Fourier Transform


What's a good classification algorithm for a large amount of classes? My requirements are these:

Fast classification of vectors of length 100 to one of 30000 classes
Iterative learning (can improve the model after it was first learnt)
Preferably available implementations in Matlab and Java

Thanks

I would suggest using Deep Learning and especially convolutional neural networks. This might be a generic approach to large scale classification tasks. One such example is identifying faces. As an example, 10.000 classes have been tackled by Sun et. al. in the following : http://ieeexplore.ieee.org/document/6909640/

And iterative learning is of course supported.


What is the official definition of a feature in signal classification? Nowadays read any paper or reference on signal classification scheme you would invariably come across the term "feature". What does feature mean exactly? What is considered a time domain feature, what is considered a frequency domain feature and what is considered a phase domain feature?

I would be glad if someone can clarify this concept for me once and for all.

Features are properties of the signal used to describe it in some way. Sometimes the term "descriptors" is used, for example MPEG-7 audio-visual description standard:

Time domain features are computed from the audio waveform, while frequency-domain or spectral features are derived after performing an FFT.

You should look through a feature-extraction software package like Yaafe, to get some practical intuition. There are also many tutorials and papers online for further learning.


What are the largest FIR filter lengths used in practice? This is a rather general set of questions.

I'm interested in knowing up to what orders, direct form FIR filters are considered useful for practical applications (both in hardware and software)?
What if fixed-point arithmetic is used?
Are FIR filters with several hundred or thousand fixed-point taps suitable for some problems, even if this is rare?

All applications are bound to some notion of accuracy of the magnitude/phase response.

Given a FIR is a polynomial approximation to a desired magnitude and phase response:

The lower the number of taps, the poorer the approximation, but the
lower the computational overhead.
The higher the number of taps, the better chance you have of a more
accurate approximation (depending on your design approach), but the
higher the computational overhead.

There's not much more can be said on the finer points unless you ask a more specific question at which point there are many people cleverer than me in this forum who will be able to give you some good pointers.

For instance FFT makes high order FIR filters possible as it is computationally efficient and has good numerical properties. High tap FIR filters aren't rare and I can vouch for them being used in audio (64K and 128K taps are common for some purposes).

Using fixed point arithmetic usually implies you are computationally bound by a chip without an efficient floating point arithmetic unit and the numerical properties of fixed point arithmetic tend to limit the accuracy of high order filters. So high order filters are not aided by fixed point arithmetic.

UPDATE, based on question in comments about upper limit on taps for fixed point:

I've only ever implemented FIR filters using floating point arithmetic, but I have written fixed point routines (for non FIR related applications) on a 32-bit 16 instruction RISC architecture in assembler. It's possible to get any precision you want by extending your arithmetic across words at the expense of clock cycles i.e. you can improve your precision in fixed point if necessary, but at some point you will hit the limit of your processing capability. So you ideally need to say what chip you're using, what language/libraries you're going to use, show us a graph of your ideal magnitude/phase response, what latency you're happy with and what input/output rate. It's impossible to give hard limits without some real constraints. It might be another question on SO perhaps...


Convert image of an irregular surface to it's signal equivalent I want to take the image of an irregular object , say

and then I want to convert the  elevation of the points into a spatial signal. To make my point clear: the signal would be 'height of the object's surface from an assumed datum' vs 'distance'.

For example

this is the kind of signal that I want to construct from my image.

What software/tool can I use to convert my image to it's signal equivalent?

As I can see from comments, you are taking multiple images from the same camera at diferent angle. So what you might be looking for is Structure from motion .This is similar to Stereo vision. You can find open source solutions and software tools available here. Good luck.


Where is the "time and frequency" in a Discrete Wavelet Transform? The DWT is motivated by the ability to see signal in time and frequency at the same time, if I understood that correctly.

I am using the following Matlab commands for some signal x:

[A,D] = dwt(x,'db4');

Then I plot A vs some vector t, representing 'time'. I do know the algorithm behind A and D: low (high) pass filter, then downsample by 2.

My questions: Where is the time axis and where is the frequency axis if I all what I have are the coefficients A (approximation) and D (detail)?

Second question: If all of this is a simple filtering and downsampling, I don't see what is new in DWT...

DWT is motivated by the ability to analyse the signal in time and frequency with an adaptative time and frequency resolution.

Some clever algorithm allows you to perform a wavelet analysis with a recursive implementation. Matlab's DWT does that, matlab's dwt perform only one level of DWT, therefore, A is only the output of the low-pass filter while D is the output of the high-pass filter. There is no time and frequency here.
It only gives you a very coarse information about the time-frequency content of the signal.
Please check how to implement a DWT.

Second question : To see the novelty, you have to compare that to short-time frequency transform. With the DWT the time-frequency resolution depend on the frequency while with the STFT it's always the same. This is a VERY DESIRABLE property : You want a short time windows for quicly mooving high frequency signals and a long window for slower low frequency signals.

Moreover, please consider that the filtering/downsampling implementation is a very refined and efficient one. DWT is an incredible tool for analysis, denoising, feature extraction etc ...


Is it possible to do further beamforming on already beamformed data? I've a circular array of eight FM antennas. I'm doing digital beamforming to form eight beams with the following beam patterns:

The objective is to reduce my side lobes even further, and to get narrower main beam. What I'm interested in to know if it is possible to take two adjacent beams, and do another step of beamforming, i.e., considering them as if coming from two antenna array with the beam pattern shown in the figure. I don't know how to incorporate the beam pattern into the standard set of equations for beamforming, as I think an omni-directional antenna is assumed. Any directions, suggestions are welcome.

For now I will focus on the conceptual, pending feedback from you:

  What I'm interested in to know if it is possible to take two adjacent
  beams, and do another step of beamforming,

The answer to this is yes, this is possible to do, and this is done in many beamforming applications. Think about it like this: All beamforming is really doing, is undoing the delays your signals exhibit (by virtue of different sensors at different locations), and summing the result. (Example: Delay-And-Sum beamformer). But in what order you undo the delays, before summing doesn't matter.

Thus, the beamformed result of a sub-group of sensors can be considered as one 'element', and the beamformed result of another sub-group of sensors, another 'element', and then both 'element's' outputs are beamformed.

In other words, first we remove the effects of the delays in one sub-array. Then we remove the effects of the delays in the other sub-array. At this point, both signal outputs might still have a delay relative to each other. So then again, both those inputs can be beamformed again. You can keep doing this ad infinitum.

  The objective is to reduce my side lobes even further, and to get
  narrower main beam.

If this is the objective, then the only way to reduce your main lobe width is to sample the spatial field with more sensors. This is the direct analogue to time-domain sampling. If we want greater frequency resolution (smaller main lobe width in the DFT analysis), we need to process a greater number of samples in time. Similarly, if we would like a smaller main-lobe width of our spatial spectra, we need to sample the spatial field with more sensors.

EDIT: Taking from the above, if you want to decrease your mainlobe width, you need more spatial samples. However if you want to reduce your sidelobes, then you would use a variety of windowing/shading functons. (Each of those functions will give you a slight increase in mainlobe width, but that is the price to pay). I am not sure what shading you have used for the above plots, but that will be good to know.


How can a signal have no DC component I am reading a nice book on DSP

  The Scientist and Engineer's Guide to Digital Signal Processing

It mentions

  In electronics, the mean is commonly called the DC (direct current) value.

How can a signal have no DC component if the mean is the DC component? Surely every signal's values can be averaged?

A signal having a mean-value or DC component of zero is commonly referred to as mean-free or as having no DC component. It does not mean that it cannot be averaged, just that the average comes out as zero.

Might be a little inexact but it is very common.


Image Contrast Enhancement Threshold Selection I have a a set of 24bit images with bright reflective spots washing out the underlying detail. While I'm interested in tone mapping algorithms, since these bright spots are purely noise I would like to threshold them out and render the rest. The histogram of the image looks like:

and is slightly bimodal in the tail where the overexposed pixels form a slight bulge. Sometimes this bulge is larger or smaller, but it's never dominate and the tail of valid pixels is fairly long.

Is there an efficient approach to identify the boundary of a bimodal distribution (even in this challenging case)? I recall an algorithm that iterates to identify a midpoint between two unbalanced "chunks" of a distribution but cannot recall its name. Ideally, the algorithm would find a theshold of about 3500 so that the very bright pixels are thresholded off.

What you are looking for called Otsu Method / Otsu Thresholding.
It solves an optimization problem to set the optimal threshold between two modal distribution.


RNN as posterior probability estimation in speech recognition with HTK I'm new to speech recognition and deep learning and in a learning phase.

I'm trying to follow this paper to learn how to use RNN as posterior probability estimation in an HTK environment. The paper proposes RNN-HMM hybrid system, so for the HMM part I need to use HTK platform.

The problem is that I couldn't even start from anywhere. I have a sample code which uses HMM to recognize digits, but I'm unable to solve at which part should I insert RNN to the code.

If there are any ideas, I would be glad.

The system in the paper is shown as follows:

I have codes in python environment and HMM is applied using HTK.
After converting data to MFCC format, I should use RNN, but after using RNN, which steps of HMM should I apply to generate RNN-HMM acoustic system.

It is very huge effort to build RNN-HMM system with HTK. It is easier to use more flexible toolkits like Kaldi, they have RNN-HMM already implemented.

HTK in latest version 3.5 beta already has support for DNN (not RNN but still a neural network), you can grep for HYBRIDHS in sources. For details on DNN-HMM it is better to read HTKBook 3.5 section 3.12 DNN-HMM Systems and Chapter 14 Artificial Neural Networks

You need to modify all places in HModel.c file to introduce your own type of the model. It is not a separate step you could introduce, you have to modify complex interaction between acoustic model and HMM search code.

It is also not straightforward to embed Python code into HTK C code, you might dump RNN scores in a file and read them from C, otherwise, you will have to reimplement your RNN code in C.


What is the computational complexity of Spectrum Estimation Algorthms and What they depend upon? There are lots of spectrum estimation techniques, each with some pros and cons.
 Algorithms like:

Music
Welch's method
Yule-Walker AR method
Periodogram
modified covariance method
multitaper method (MTM)
Spectrogram using short-time
Fourier transform
Burg method covariance method

and many many other.
 However I do not know which is more proper for embedded system implementation, because I don't know their computational complixy so I am more interested in mathematical estimation of computational complexity of this methods and especially how their complexity is affected with factors like input noise, spectrum occupation, etc?

I'v found Spectrum Estimaton Methods, which greatly explains these methods, but does not compare them from computational complexity point of view.

The computational requirements for various non-parametric PSD estimation methods is discussed in [*]. The complexity depends on the length $N$ of your data record. Roughly speaking:

The most straightforward periodogram method which calculates the magnitude squared of the DFT of $N$ samples has a complexity of $N \log N$.
Bartlett and Welch methods use $M$ overlapping windows of length $M$ each and calculate ~$N/M$ DFTs, so the complexity is $\frac{N}{M} \cdot M \log M = N \log M$.
Yule-Walker method for an $AR(p)$ inverts a $p\times p$ Toeplitz matrix which has a complexity of $O(p^2)$. The Toeplitz matrix itself is formed from estimates of the autocorrelation values calculated at $p$ different lags from the $N$-long data vector giving a complexity $O(Np)$. So the total complexity is $O(p^2 + Np)$.
MUSIC estimates an $M \times M$ autocovariance matrix from data by averaging the outer products of, say, $N/M$ $M$-long blocks of the original $N$-long signal and then calculates the eigendecomposition of the $M\times M$ autocovariance matrix which is $O(M^3)$. ($M$ is chosen by the user. For a signal with $p$ frequency components, choose $M>p$.)

I don't think the complexity is affected by input noise level because the algorithms themselves do not change.

Reference

[*] J. Proakis Digital Signal Processing, 4th Edition Ch. 14
https://www.pearsonhighered.com/program/Proakis-Digital-Signal-Processing-4th-Edition/PGM258227.html


Can the DFT can be performed without normalization? If we have OFDM system with $N$ sub-carriers, the DFT matrix can be expressed as follows:
$F = dfmtx(N)/sqrt(N)$
My question, is it possible and practical to use the matrix without Normalization, it means use the iDFT as $F = (dfmtx(N))'$, and then at the receiver end, I will multiply with  that matrix and divided by $N$  (like this $dfmtx(N)/N)$ to get scalar diagonal matrix whose diagonal elements are 1?
thank you

Yes.  In fact, depending on your receiver architecture, you may not need the divide by $N$ in the receiver.
Generally, in a receiver you care about the relative phases of the bins (or possibly their strengths with respect to some other part of the same signal).  So there's no reason to do that multiply if you don't want to.
If you're working with fixed-point math, just make sure that you don't experience over- or under-flow.


Suggestions for FFT based tools to analyze music - Preferably MATLAB/free software packages I am interested in analyzing full length songs for general features of the track which are not otherwise apparent. For e.g., I want to look at various parts of the song with time-frequency plots/wavelets/other signal processing tools for, lets say -  How does the frequency content of the song during the guitar solo compare to the part when only vocals are on? What is the frequency with which the guitar solo starts? What is its peak? ..and so forth.

Although I use FFT often for research, this is something I am fascinated by.... Any suggestions where to start?

At this point I have got WAV files for MATLAB analysis, since they have almost no loss. Where do I go from here? Any references to concepts/open source codes/literature etc. would be greatly appreciated.

Also, I am not into audio engineering or anything similar. Its just for my own interest.

My favorite is Sonic Visualiser. It is free and fast. In combination with the Vamp Plugins you can analyze low level features up to high level features like onsets, pitches, chords and so on. It has also an nice spectrum and spectrogram fea


Covariance matrix associated with random DC level in Gaussian noise Given a signal

$x[n] = A + w[n]$

where $A$ is a Gaussian random variable and $w[n]$ is Gaussian white noise, then the covariance matrix of the signal is given by

$[C(\sigma^2_A)]_{ij}=E[x[i-1]x[j-1]]=E[(A+w[i-1])(A+w[j-1])]$.

I understand this definition, but where does the reduced expression

$=\sigma^2_A+\sigma^2\delta_{ij}$

come from?

It is important to add that $w[n]$ is uncorrelated with $A$. This means that

$$E\{Aw[n]\}=E\{A\}E\{w[n]\}=0\tag{1}$$

because we may assume that $w[n]$ is zero mean: $E\{w[n]\}=0$. In this case all cross-terms in the expectation are zero, and we're left with

$$C_{ij}=E\{A^2\}+E\{w[i-1]w[j-1]\}\tag{2}$$

Since $w[n]$ is white, you have $E\{w[i-1]w[j-1]\}=\sigma^2\delta_{ij}$. The expectation of $A^2$ is given by

$$E\{A^2\}=\sigma_A^2+\mu_A^2\tag{3}$$

where $\mu_A=E\{A\}$ is the mean of $A$.

So we finally have

$$C_{ij}=\sigma_A^2+\mu_A^2+\sigma^2\delta_{ij}\tag{4}$$

The expression in your question is only valid if $\mu_A=E\{A\}=0$, i.e., if $A$ is a zero mean random variable.


Online DFT Algorithm I have a discrete audio stream $x$ that needs to be processed in real-time. Specifically, as the each new sample is received, I would like to compute a Fourier transform of the last $n$ samples of the signal. So, for example, upon receiving the $i^{th}$ sample, I would like to compute the Fourier transform of the sequence $(x_{i-n+1}, x_{i-n+2}, \dots, x_{i-1}, x_{i})$. However this means that when I receive the $(i+1)^{th}$ sample, I will be computing the Fourier transform of $(x_{i-n+2}, x_{i-n+3}, \dots, x_{i}, x_{i+1})$, which is almost identical to the first sequence. So my question is, having already computed the Fourier transform of the first sequence, is there some way to compute the new Fourier transform from that, instead of having to recompute it from scratch each time?

(Note: the paper pointed by hotpaw2's link is actually describing in more detail the algorthm I presented here)

Consider a data window length of $N$ samples from $n=0$ to $n=N-1$. Let your original data window be $x_1[n]$, whose first sample is $x_{old} = x_1[0]$. Now your new data set is denoted as $x_2[n]$ whose samples are actually one sample left shifted verisons of $x_1[n]$, i.e. $x_2[n] = x_1[n+1]$ for $n=0,1,...,N-2$ plus a new arrived data to position $n=N-1$ which is denoted as $x_{new} = x_2[N-1]$.

Then the following algorithm will compute the N-point DFT, $X_2[k]$ of the new data set $x_2[n]$ from that of the already computed and stored N-point DFT $X_1[k]$ of the old data set $x_1[n]$ as:

$$ X_2[k] = e^{j \frac{2\pi}{N}k }( X_1[k] + (x_{new} - x_{old}) ) $$

for each $k=0,1,...,N-1$.  This updated computation of $X_2[k]$ from pre-computed $X_1[k]$ requires N complex multiplications and N real additions. Compared to a direct N-point FFT which requires $N \log_2(N)$ complex MACs, this would therefore be an improvement by a factor of $\log_2(N)$ which for example would render to a factor of 10 when $N=1024$ using low level languages such as C.


FFT (DFT) of 2D Images Branching off from the question [Mapping the value of a sample in a 2-D DFT to cycles/pixel], I have a few more questions related to 2D FFT of images.

I correct in assuming that the fundamental frequency is the image itself?
And one copy of the image is equivalent to 1 wave cycle?
The sampling frequency is the number of pixels in the image?
And If this is correct, how do I compare two distinct images?
Can I compare the sampling frequencies ( $\frac{1}{N_1}$ & $\frac{1}{N_2}$ ) and frequency peaks of one image with a second image ($\frac{1}{N_2}$)? ($\frac{1}{N_1}$ & $\frac{1}{N_2}$ are the dimensions of an image)

What I am trying to understand is - what the X & Y axis of a complex 2D FFT plot is? and how to compare information from 2 plots?

Consider a 1D wave of 10sec duration.

I correct in assuming that the fundamental frequency is the image itself?
Does the 10 sec sample equal the fundamental freq? No. same for images.

And one copy of the image is equivalent to 1 wave cycle?
Repeat above. Answer: No

The sampling frequency is the number of pixels in the image?
One pixel can be thought of as the sampler. But, this could be changed. You could go for sub-pixel sampling also.

With this I think the last question has no relevance. Comparing two images, to the best of my knowledge is an open ended problem.


DFT of the same signal with different values of N Let $x[n]$ be a discrete signal of 2 samples. We know that its DFT with N=4 is $X[k]=[0, 1+j, 2, 1-j]$. Without calculating $x[n]$, how can we know the DFT with N=2?

I have tried to use the relation between DFT and DTFT. Specifically, that the DFT is the result of sampling the DTFT in the frequency domain every $\frac{2\pi k}{N}$ samples. However, I don't know how to continue. I would really appreciate some help.

$$X_4[0] = x[0] + x[1]$$
$$X_4[2] = x[0] - x[1]$$

Now $$X_2[0] = x[0] + x[1]$$
$$X_2[1] = x[0] - x[1]$$

Therefore $$X_2[0] = X_4[0]$$
 $$X_2[1] = X_4[2]$$

An alternative explanation is the following, if we were to decimate in the frequency domain by 2 (the 4 point FFT), you would incur circularly shifted Aliasing in discrete time domain. Since we started with a two point sequence then zero padded it to get a 4 point sequence, this decimation will result in Aliasing at those zero stuffed places leaving your original sequence non aliased. So say if we were to decimate more than 2 the sequence would incur Aliasing in time domain. So the fact that we are still able to have the original 2 point sequence non aliased after the decimation, therefore the FFT values of this 4 point FFT still carries the information about the 2 point FFT

Ex: if say we would have a  4 point sequence and you zero pad to get a 16 point sequence and take the FFT, now if you decimate the FFT by a factor of 2 you would get Aliasing in samples 8 to 15 of the time domain sequence.


What are the mathematical consequences of a discrete signal living in Hilbert space? I am doing the Coursera course on DSP by Prandoni and Vetterli, and was excited to learn that we can use the tools of vector spaces to analyze discrete signals. At this early point in the course the professors did mention that since the signals live in Hilbert space, we can naturally talk of orthogonality, for one.

I am curious: what other linear algebra operations/concepts do we use with discrete signals, and what is the rationale behind them?

whether or not a signal lives in a Hilbert space depends only on if you say it does.  if you say so, then you need to specify an inner product $<x,y>$ (which is pretty easy for finite-length discrete-time signals)

$$ \langle x,y \rangle = \sum_n x[n] y^*[n] $$

and you need to use that inner product definition to define a norm

$$ |x| = \sqrt{\langle x,x \rangle} $$

which should also be the same as the distance measure from $x$ to the $0$ element in that Hilbert space.

and you need to be able to add and subtract elements in the Hilbert space (easy to do for discrete-time signals).  why this might be useful would be for defining a signal $x$, an estimation of that signal $\hat{x}$, and an error signal, $e_x$ such that

$$ \hat{x} = x + e_x $$
or
$$ e_x = \hat{x} - x $$

we can think of the norm of $e_x$ which is

$$ |e_x| = \sqrt{\langle e_x,e_x \rangle} $$

as a measure of the distance that $\hat{x}$ is from $x$.

if gives you a formal framework for defining error and then maybe help you with algorithms (like a matched filter) for minimizing error, so that a received signal most closely approximates what was transmitted.  or so that some other  "thing" most closely hits some specified target "thing".


How to filter a signal using a bandpass filter consisted of two moving average filters? I want to filter a PPG signal on a microcontroller. I have limited memory and a was searching for low computational methods. I found the work of Kazuhiro Taniguchi, Earable POCER: Development of a Point-of-Care Ear Sensor for Respiratory Rate Measurement where they use moving average filters (m3, m30 and m80) in order to filter their data and obtain freequencies between certain values, specifically they create a a passband between 189 mHz and 504 mHz. Their idea is that after the initial m3 moving average(every 3 values) they create the m30 and m80 and they apply an iteration between them ( r= m30 - m80). This is their way to obtain the values corresponding to the passband of interest, eventhough, as they specify,

The moving average does not have an ideal lowpass filter
function, and so some frequency elements other than those in the passband may pass, even if they are attenuated.

to obtain the window size for their moving average filters they applied the following equation:

I couldn't quite understand where this equation (5) was coming from but as soon as you replace the value $n$ with 30 or 80 you get the values of the passband cut-off freequencies 189 mHz and 504 mHz.
When I asked for specifications on how they did it, they send me to a japanese forum that unfortunately I couldn't translate but there were two links to stackexchange, on cuttoff freequency and filter design.
I tried to adapt all these new staff on my model but couldn't get the passband wanted (0.1 Hz - 0.8 Hz) using the equation (5) from the image above with my parameters (sampling freequency of 50Hz, cut-off freqs etc...).
I don't understand what is the problem and my question is the following:
What are the moving average filter windows that I have to use and in which order to be able to filter a signal with a sampling frequency of 50 Hz in order to isolate the frequencies between the passband 0.1 Hz - 0.8 Hz?

The Japanese link actually implies how to derive Eq.5, and ironically they also refer to an existing dsp.se answer at the bottom.
Derivation of Eq.5 is as follows:
Consider a moving average filter of length $N$, with the impulse response $h[n]$:
$$h[n] = \begin{cases} ~~~1/N~~~,~~~n=0,1,...,N-1 \\ ~~~0~~~,~~~ \text{otherwise} \end{cases} \tag{1}$$
Magnitude of its frequency response $~H(\omega)~$  (DTFT of $h[n]$) can be shown to be:
$$ |H(\omega)| = \frac{1}{N} \left|\frac{ \sin(\frac{\omega}{2}N)}{\sin(\frac{\omega}{2})}\right| \tag{2}$$
Now, after replacing the $\sin()$ functions with their Taylor expansions in powers of $\omega$, it can be shown by polynomial long division that Eq.2 is also given by :
$$|H(\omega)| = 1 + \frac{1}{24}(1-N^2) \omega^2 + H.O.T. \tag{3}$$ where H.O.T. refers to higher order terms in powers of $\omega$. An approximation for small values of $\omega$ is obtained by neglecting H.O.T.:
$$|H(\omega)| \approx 1 + \frac{1}{24}(1-N^2) \omega^2 \tag{4}$$
Using this approximate frequency response magnitude, we can obtain an approximate cutoff frequency $\omega_c$ at which the magnitude $|H(\omega_c)|$ falls to $1/\sqrt{2}$ of its value at $\omega = 0$, which is $H(0) = 1$.
$$ |H(\omega_c)| = \frac{1}{\sqrt{2}} = 1 + \frac{1}{24}(1-N^2) \omega_c^2 \tag{5}$$
Replace the discrete-time frequency $w_c$ by $w_c = 2\pi f_c /f_s$, where $f_c$ is the analog cutoff frequency in Hz, and $f_s$ is the sampling frequency in Hz. Finally solving the resulting algebraic expression for $f_c$ yields the formula that you refer to as Eq.5 in the document :
$$ |H(\omega_c)| = \frac{1}{\sqrt{2}} = 1 + \frac{1}{24}(1-N^2) \left( 2\pi \frac{f_c}{f_s} \right)^2 \tag{6}\\\\$$
$$ f_c = \frac{1}{\pi} \frac{\sqrt{6 - 3\sqrt{2}}}{\sqrt{N^2-1}} ~f_s ~~ =~~  \frac{0.422}{\sqrt{N^2-1}} ~f_s \tag{7} $$
Eq.7 above is the formula that provides an approximate cutoff frequency calculation for the moving average filter of length $N$ (order $N-1$). In your posted link, Eq.5 there's a slight variation at the scale of $0.442$ instead of $0.422$, probably they tried to apply some correction to the actual cutoff vs approximated one.
Note that in the derivation, we've used an approximation of the DTFT magnitude which was valid as long as $\omega$ was small compared to $\pi$. This means that the approximation will be satisfactory if $\omega_c$ is close to $0$, or in other words, $f_c$ is a small compared to $f_s$. And indeed this will be the case for high order moving average filters. And the approximation gets better as $N$ increases.
Using such two moving average filters with approximate cutoff frequencies $f_{c1}$ and $f_{c2}$ to create a bandpass filter will not be very satisfactory unless your out of band signal energy is insignificand after 20 to 30 dB of attenuation.


The cost function of LS algorithm I'm trying to understand the derivation of LS algorithm for channel estimation in OFDM.
Could anyone explain why do we need complex conj transpose in the below equation?

Because when a quantity can be complex, and even when it is just real, the absolute squared difference $|f-g|^2$ can be expressed in both domain (complex and real) as:
$$|f-g|^2 = (f-g)^H(f-g)$$
and of course this is correct as well for reals. This setting is often related to Hilbert spaces.


Audio EQ Cookbook without frequency warping The famous http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt offers a set of [real] biquad filter calculation formulas that generally work fine.

However, when filter's frequency approaches Nyquist frequency the Q (bandwidth) specification of a filter becomes greatly distorted - usually it shrinks a lot (even though the author mentioned that he performed a necessary pre-warping).

I'm in quest for filter formulas that do not have such strong bandwidth distortion. I need peaking/bell, band-pass, low-pass, high-pass, high-shelf, notch filters. I know this can be done as previously I've bought a peaking/bell/band-pass filter formulas with less distortion, but they are still not perfect and I need other filter types.

So, I'm also willing to pay for the solution if the price is right.

Alternatively, if one could point me to an optimization algorithm that works with Z-domain filters that would be great, too. Unfortunately, most usual optimization algorithms do not work well in Z-domain - they can't optimize a set of parameters to match a desired frequency response (probably due to periodic functions used to calculate the frequency response).

here is a quick look at how the 5 degrees of freedom for the parametric EQ can be viewed.  it's my take on what Knud Christensen of tc electronic came up with about a decade ago at an AES convention and this patent.

so, forget about the Cookbook (and the issues of Q and bandwidth therein) and consider (in the s-plane) the parametric EQ as the sum of a bandpass filter (with a $Q$ value) in parallel with a wire:

$$ H(s) = (G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + 1 $$

$G_\text{boost} = 10^{\frac{dB}{20}}$ is the gain of the peak (or valley, if $dB<0$).  the gain at DC and at Nyquist is 0 dB. that's a 2nd-order IIR and there are 3 independent parameters. 2 more to go.  so we next add an overall gain parameter:

$$ H(s) = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + 1 \right) $$

that's 4 knobs to twist.  one more knob to add (without raising the filter order) and we will be done with adding more independent parameters.

so what Knud does here is replace that "wire" (that trailing "$1$" in the transfer function) with a prototype shelving filter that must have the same poles, the same $Q$ and $\omega_0$ as the BPF, so that the denominator is the same.  the transfer function of that shelf is:

$$ H_\text{shelf}(s) = \frac{R\left(\frac{s}{\omega_0}\right)^2 + \frac{\sqrt{R}}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} $$

where $R \triangleq 10^{\frac{tilt}{20}}$ and $tilt$ is the gain differential of the shelf in dB.  this is what offsets the gain at Nyquist to be different than the gain at DC.  After bilinear transformation, Nyquist gets boosted by $tilt$ dB and the gain at DC remains unchanged.  Like the $dB$ boost parameter, the $tilt$ parameter can be either positive or negative.  $G_0 \cdot R$ is the linear gain at Nyquist.

put that all together and you get:

$$ \begin{align}
H(s) & = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + H_\text{shelf}(s) \right) \\
 & = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + \frac{R\left(\frac{s}{\omega_0}\right)^2 + \frac{\sqrt{R}}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} \right) \\
& = G_0 \frac{R\left(\frac{s}{\omega_0}\right)^2 + (G_\text{boost}+\sqrt{R}-1)\frac{1}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1}  \\
\end{align}$$

no matter how you look at that, this has 5 degrees of freedom and those 5 biquad coefficients are fully defined from these 5 parameters.  doesn't matter if you map from $s$ to $z$ using the blinear transform or the trapezoidal rule (effectively the same thing) or any other method that does not change the order of the filter.  you may have to fudge the definition of $Q$ or bandwidth, you may have to compensate $\omega_0$ and/or $Q$ for frequency warping effects (like you get with the blinear transform), but if you paid big bucks for something that gets you a 2nd-order IIR filter, it doesn't matter if you implement it with any Direct Form or transposed Direct Form or Lattice or Normalized Ladder or Hal Chamberlin's State-Variable or Andrew Simpson's modeling of linear analog with trapezoidal integration, eventually you get to 5 coefficients and they can be mapped to these 5 independent parameters.  it's all the same.  whether or not you paid money for a license or not.  the math is stronger than any claims made by whomever you are licensing from.

Just FYI, i solved where the true peak or valley frequency is when there is a $tilt$ that is not zero.  the frequency where the peak or valley has been nudged over by the tilt is:

$$ \omega_\text{peak} \ = \ \omega_0 \
\sqrt{ \frac{Q^2 \left(R - \frac{1}{R}\right)}{G_\text{boost}^2 - R + 2 Q^2 (R - 1)} + \\ \sqrt{\frac{G_\text{boost}^2 - \frac{1}{R} + 2 Q^2 \left(\frac{1}{R} - 1\right)}{G_\text{boost}^2 - R + 2 Q^2 (R - 1)} + \frac{Q^4 \left(R - \frac{1}{R}\right)^2}{\left(G_\text{boost}^2 - R + 2 Q^2 (R - 1)\right)^2} } } $$

you can see that when $tilt = 0$, then $R=1$ and consequently $\omega_\text{peak} = \omega_0$.   the peak gain $G_\text{boost}$ might also have to be adjusted a little and that has yet to be worked out.  a good first guess would be $G_\text{boost} \leftarrow \frac{G_\text{boost}}{\sqrt{R}}$ or maybe $G_\text{boost} \leftarrow G_\text{boost}-(\sqrt{R}-1)$.


which spatio-temporal properties does dynamic textures exhibits? know that dynamic textures are spatially repetitive and temporally varying patterns.

1.I read one sentence in one paper i.e. They are sequences of frames exhibiting certain stationarity properties over time. They exhibit spatial self-similarity and temporal self-similarity. but i don't understand how it does?

2.Also how examples like "traffic on road" or "rafting" can be considered as Dynamic Textures? Because i don't thing they exhibit any stationarity property over time.is it?

Stationarity over time just means that if you come back and look at the traffic sometime later it will look pretty much the same as it did when you first saw it. Sure there will be different cars but even the new cars will look similar to the first ones. The cars might be in different positions but their random arrangement will have the same properties as the original arrangement.


'imread' function in MATLAB gives image data.what this image data is?is that intensity of pixels?and how MATLAB computes this data? i checked imread function code but couldnt understand how it is calculating the data from given image.Please explain.

'imread' function just read the RGB(in case of color image) or black/white(in case of black and white image) value of each pixel in image and shows it in a matrix format.
it is not 'calculating' these values its just reading them from the image/path you provide to imread function.


Error plot between known and estimated data The channel is an FIR model with input $u$. The input takes in values which are symbols from some constellation. Using an equalizer such as the Least Mean Squares (LMS), I estimate the input to the channel from the known noisy observations and the channel coefficients. Let the estimated input to the channel be denoted by $\hat{u}[n]$.

The coefficients of the channel are denoted by $c$ of length $L$. The noisy output of the channel is $y$. Therefore we have, $y[n] = \mathbf{c}^H*\mathbf{u} + v[n] \tag{1}$ for $n = 1,2,\ldots,N$ where $v[n] \sim N(0,\sigma^2)$ is a Gaussian noise. The weights of the equalizer are given by $w$ and the equalizer is of length $2L+1$. The input to the equalizer is $y$. The output of the equalizer is given by $\hat{u}[n] = y*\mathbf{w}$ where $*$ denotes the convolution. My questions are:

1) Mean Square error between the estimated input and known / desired input to the channel for different values of Signal-To-Noise ratio (SNR) --The result of equalization will give us the estimates of the input to the channel. The estimate of the input to the channel is $\hat{u}[n]$. If I need to calculate the error between the estimated input and the known input, I should probably use the formula: $MSE = \frac{1}{T}\frac{\sum_{n=1}^{N}\hat{u}[n] - {u}[n]}{N}$ and $T$ denotes the number of independent trails. As an example, the input can take values from the set $\{1,2,3,4\}$. The first 10 symbols of the data can be $2,1,3,1,1,2,4,1,2,3$. These are the known values. Then to validate the model after equalization,let the estimates input symbols to the FIR channel be $\hat{u} = 1,3,1,1,1,2,4,1,2,2$. The error, $e = (u[1]-\hat{u}[1])^2+(u[2]-\hat{u}[2])^2+(u[3]-\hat{u}[3])^2 + (u[4]-\hat{u}[4])^2 +(u[5]-\hat{u}[5])^2 + ...+(u[10]-\hat{u}[10])^2$
Is this the correct way?
I am not using the communication tool-box in Matlab, and building my own code

2) Symbol error rate -- In communications, I have seen the symbol error plot (please refer :http://www.dsplog.com/2008/03/29/comparing-16psk-vs-16qam-for-symbol-error-rate/).

How is Using the same data which is used to obtain the MSE plot between the estimated symbols and known symbols, how can I obtain the SER plot? An example in Matlab will be very helpful.

Thank you.

Your formula/method for computing MSE between estimated and known inputs looks good to me. For symbol error rate you could use something like a Hamming distance which simply counts the number of times the estimated symbol is different from the actual symbol. In your 10 symbol example the error rate is 4/10 i.e 40%.

In Matlab you can do something like:

symbol_error_rate_percent = sum(u_hat ~= u)/length(u) * 100

You can compute this for $T$ trials and average them, like you did for MSE. Or you could just make a really long sequence of length $10 \times T$ and do it once.


Representing stereo signal with complex numbers as input to DFT Usually, when we have a stereo signal, we process each channel separately in order to extract the frequencies using Fourier transform.

However, Fourier transform can also be applied to complex numbers as well. So what if I represent the original stereo signal with complex numbers, do a Fourier transform, and back?

I know that in this case, for DTF, I would need to keep all N frequency numbers as opposed to N/2 for DFT on the real data. Should I then be able to recover my original signal from the frequencies? (My guess would be "yes".)

Now, what if I want to detect if a specific frequency is present in the signal? I could use the formula for DFT:

So, what if instead of real $x_n$, I pass a complex number containing values of the stereo channels as real and imaginary parts? Would that make any sense? What would that detect?

For your first claim. Yes it's true. You should keep the full N-DFT complex frequency samples as opposed to half, N/2, for the real input case. And also yes, you will retain back the original stereo signal perfectly by an inverse DFT.

Whether this artificially produced complex signal adds any new information that cannot be obtained from individual channels is not very clear, however. At least from the linearity of DFT we can guess that the resultant DFT of the complexed-stereo pair will be exactly given by those individual left/right channel DFTs, and this suggests to me that, you wont be obtaining some interesting benefit here.

Keep in mind that, the (geometric) spatial relationship between a pair of stereo microphones can be investigated by a joint processing of those channels. However, a DFT alone does not care for any spatial relationship between real and imaginary part of its time domain complex input. Hence I suspect any good use of that for spatial processing.

Yet, in some cases, there may be a computational benefit due to performing a single complex input DFT instead of two real input DFTS. In such cases you are favored to do so.


Second Order Transfer Function : Imaginary component in $R$? I have been looking at the following transfer function:

$$ H(z) = \tfrac12 -  \tfrac12 z^{-2} $$

Given the usual method for finding $\theta$ and $R$ in the complex plane, I calculate that $\theta = \frac{\pi}{2}$ and $R = j\sqrt{\frac12}$

This would initially make me think that the conjugate zero pair here is 90° from the real axis, however this imaginary $R$ component has thrown me completely - the equivalent frequency response also seems odd to me as it shows zero at DC and Nyquist with a gentle curve upwards in between.

What is going on here? I think that $R$ in the complex plane can't be imaginary, but don't understand how to arrive at the correct angle and magnitude for the conjugate poles.

HINT:

If you're looking for the zeros of $H(z)$, try $z_0=1$ or $z_0=-1$. If that works, try to figure out why, and derive it yourself. Then, think about what the magnitude and the phase of those zeros might be.

EDIT (in reaction to your comment): If you have a second-order polynomial with real-valued coefficients, then two complex conjugated zeros are not the only option; you can also have two real-valued zeros.


Running window maximum: Any cautions? Are there any statistical concerns with using a running-window maximum to capture peak values?

Background: We are measuring acceleration during random impacts on a structure. The impacts need ~500Hz sampling to resolve, but due to data volume I can only permanently save a 10Hz data vector. One workaround is to record a running window maximum since I am most interested in peak accelerations during impacts, and how often impacts occur. Are there any statistical concerns (bias, accuracy) about using a running-window maximum?

There's a few possible problems with taking a running maximum of all the accelerometer samples within your window.

You'll miss large negative excursions. This may not be a problem; the absolute value of excursions might be enough for your application.
You'll miss two large excursions in the one time window.
The actual times the peaks occur at will be lost, or at least be inaccurate. For example, if adjacent windows have peaks at the beginning and then at the  end,  the actual time between them will be nearly two window lengths, but the reported time will be one window length
All values will tend to be positive, which will bias the "background" acceleration.

I've written some R code to simulate discarding 49 out of every 50 values, and only keeping the maximum one in a window. The raw data are the black circles in the plot. The maximum values are the green circles; the minimum values are the red circles.

R Code Below

#27212

N <- 1000
acc <- rnorm(N)/5

Nspikes <- 10
spike_indices <- runif(Nspikes,1,N)
spike_values <- runif(Nspikes,-10,10)

acc[spike_indices] <- spike_values

decimationFactor <- 50
idx <- 1
mx <- rep(0,length(seq(1,length(acc),decimationFactor)))
mn <- rep(0,length(seq(1,length(acc),decimationFactor)))

for (k in seq(1,length(acc),decimationFactor))
{
  idxs = seq(max(1,k), min(length(acc), k+decimationFactor-1),1)
  mx[idx] = max(acc[idxs])
  mn[idx] = min(acc[idxs])
  idx <- idx + 1
}

par(mfrow=c(1,1))
plot(acc)
points(seq(1,length(acc),decimationFactor),mx,lwd=5,col="green")
points(seq(1,length(acc),decimationFactor),mn,lwd=5,col="red")


The frequency response function (FRF) fails to detect the antiresonance of a system I am trying to identify a vibrational systems by computing the frequency response function (FRF) of the system when a chirp signal is applied to its input. After comparing the FRF computed and the bode plot of the transfer function I found that the resonance has been well detected but the anti-resonance has not.

The system has the following form:

$$\textrm{sys} =  \frac{s^2+\omega_{n_{z_1}}^2}{s^2+\omega_{n_{p_1}}^2}, \quad \textrm{where } \omega_{n_z} = 10\textrm{ Hz} \textrm{ and } \omega_{n_p} = 100\textrm{ Hz}.\\$$

I have a hypothesis:

The steady state response of a linear system to a sinusoidal signal is another sinusoidal signal with the same frequency with an amplitude and a phase that depend on the system. When I applied a chirp signal, the system never reached the steady state. So while the output of the system to $\sin(\omega_{n_{z_1}}\cdot 2\pi)$ should be zero, I am still obtaining some sinusoidal signal of frequency $100\textrm{ Hz}$ as a consequence of stimulating (weakly) that part of the spectrum with the transient response.

In my opinion the correct computation should be measuring only the components of the outputs that correspond to the frequency of the input at certain particular moment.

Do you know how to correct this problem?
Is it possible to apply a moving filter or window to correct this problem?

Here is my code:

% Identification of the Frequency Response of a Transfer Function with
% Resonance and Antiresonance peaks
%% Cleaning the house
clc;
clear;
close all;
%% Definition of the System
wn_z1 = 10*2*pi;           % Anti-Resonance Natural frequency 10 Hz
wn_p1 = 80*2*pi;           % Resonance Natural frequency 80 Hz
chi_z = 0;                 % Damping factor of the zeros
chi_p = 0;                % Damping factor of the poles
s     = tf('s');
sys_1 = (s^2+2*chi_z*wn_z1*s+ wn_z1^2)/(s^2+2*chi_p*wn_p1*s+wn_p1^2);
zpk(sys_1)
%% Definition of the Chirp Input
h                    = 0.001;
time                 = 0:0.001:60;
freq_ini             = 0.1*2*pi; % Initial frequency considered (0.1 Hz)
freq_final           = 240*2*pi; % Final frequency considered (240 Hz)
chirp_input          = 1.5*chirp(time,freq_ini,time(end),freq_final,'logarithmic');

frequency_evolution  = logspace(log10(freq_ini),log10(freq_final),length(time));
figure(1);
subplot(2,1,1);plot(time,chirp_input); xlabel('time (s)'); ylabel('Chirp    amplitude')
subplot(2,1,2);plot(time,frequency_evolution/(2*pi)); xlabel('time (s)'); ylabel('Frequency (Hz)');
%% Execution of the simulation
[sys_output,~]       = lsim(sys_1,chirp_input,time);
%% Computation of the FRF
fft_input   = fft(chirp_input,2^12);                                           %% Considering zero padding
fft_output  = fft(sys_output ,2^12);                                           %% Considering zero padding
fft_input   = fft_input(1:length(fft_input)/2);
fft_output  = fft_output(1:length(fft_output)/2);
fft_output  = fft_output(:);
fft_input   = fft_input(:);
FRF              = abs(fft_output)./abs(fft_input);
frequency_vector = linspace(0,1/(2*h),length(FRF));

%% Comparing the Bode Plot (Matlab) with the FRF (own implementation)
opts = bodeoptions('cstprefs');
opts.FreqUnits = 'Hz';
close all;
bodemag(sys_1,opts,'r'); hold on;
semilogx(frequency_vector, 20*log10(FRF ),'b')
xlim([0.1,frequency_vector(end)])
legend('Matlab Bode','FRF Identification')

I agree with A_A that the chirp is too fast.

I adjusted a few things in the code

increase time to $600$ (seconds)
start chirp at $0.5\textrm{ Hz}$
end chirp at $2\textrm{ Hz}$
removed zero padding (not required here)

And the result is...

The full code is

% Identification of the Frequency Response of a Transfer Function with
% Resonance and Antiresonance peaks
%% Cleaning the house
clc;
clear;
close all;
%% Definition of the System
wn_z1 = 10*2*pi;           % Anti-Resonance Natural frequency 10 Hz
wn_p1 = 80*2*pi;           % Resonance Natural frequency 80 Hz
chi_z = 0;                 % Damping factor of the zeros
chi_p = 0;                % Damping factor of the poles
s     = tf('s');
sys_1 = (s^2+2*chi_z*wn_z1*s+ wn_z1^2)/(s^2+2*chi_p*wn_p1*s+wn_p1^2);
zpk(sys_1)
%% Definition of the Chirp Input
h                    = 0.001;
time                 = 0:0.001:600;
freq_ini             = 0.5*2*pi; % Initial frequency considered (0.1 Hz)
freq_final           = 2*2*pi; % Final frequency considered (240 Hz)
chirp_input          = 1.5*chirp(time,freq_ini,time(end),freq_final,'logarithmic');

frequency_evolution  = logspace(log10(freq_ini),log10(freq_final),length(time));
figure(1);
subplot(2,1,1);plot(time,chirp_input); xlabel('time (s)'); ylabel('Chirp    amplitude')
subplot(2,1,2);plot(time,frequency_evolution/(2*pi)); xlabel('time (s)'); ylabel('Frequency (Hz)');
%% Execution of the simulation
[sys_output,~]       = lsim(sys_1,chirp_input,time);
%% Computation of the FRF
fft_input   = fft(chirp_input);                                           %% Considering zero padding
fft_output  = fft(sys_output );                                           %% Considering zero padding
fft_input   = fft_input(1:length(fft_input)/2);
fft_output  = fft_output(1:length(fft_output)/2);
fft_output  = fft_output(:);
fft_input   = fft_input(:);
FRF              = abs(fft_output)./abs(fft_input);
frequency_vector = linspace(0,1/(2*h),length(FRF));

%% Comparing the Bode Plot (Matlab) with the FRF (own implementation)
opts = bodeoptions('cstprefs');
opts.FreqUnits = 'Hz';
close all;
bodemag(sys_1,opts,'r'); hold on;
semilogx(frequency_vector, 20*log10(FRF ),'b')
xlim([0.1,frequency_vector(end)])
legend('Matlab Bode','FRF Identification')


Why does this transfer function has a second zero I'm learning about $\mathcal Z$-transforms in DSP and I have a transfer function of the following form:

$$H(z)=\frac{2-3z^{-1}}{1-1.6z^{-1}+0.8z^{-2}}$$

When I calculate zeros and poles of this system by hand,
I get these poles from the equation

\begin{align}
\frac{z^2-1.6z+0.8}{z^2} &= 0\\
p_1&=0.8-0.4i\\
p_2&=0.8+0.4i
\end{align}

And a single zero from the equation

\begin{align}
\frac{2z-3}{z}&=0\\
z_1&=1.5+0i
\end{align}

However, when I use Wolfram Alpha to compute the poles and zeros, it also lists a second zero positioned at origin. Or in MATLAB:

The question is, where does that second zero come from?


Help understanding Hough transform I'm trying to get a Hough transform to work in MATLAB, but I'm having problems. I have a really bad way of detecting peaks that needs to be fixed, but before that I need to be able to reverse the hough transform to create the lines again properly. This is the type of stuff I'm getting right now:

looks like its rotated by 90 degrees, but I'm not sure why. I'm not sure if it's my Hough space that's wrong, or if it's the way I de-Hough and draw the lines. Also could someone help improve my peak detection? The images used in the code are here

%% load a sample image; convert to grayscale; convert to binary

%create 'x' image (works well)
a = eye(255);
b = flipud(eye(255));
x = a + b;
x(128,128) = 1;

%image = rgb2gray(imread('up.png')) < 255;
%image = rgb2gray(imread('hexagon.png')) < 255;
%image = rgb2gray(imread('traingle.png')) < 255;
%%% these work
%image = x;
%image = a;
image = b;

%% set up variables for hough transform
theta_sample_frequency = 0.01;
[x, y] = size(image);
rho_limit = norm([x y]);
rho = (-rho_limit:1:rho_limit);
theta = (0:theta_sample_frequency:pi);
num_thetas = numel(theta);
num_rhos = numel(rho);
hough_space = zeros(num_rhos, num_thetas);

%% perform hough transform
for xi = 1:x
    for yj = 1:y
        if image(xi, yj) == 1
            for theta_index = 1:num_thetas
                th = theta(theta_index);
                r  = xi * cos(th) + yj * sin(th);
                rho_index = round(r + num_rhos/2);
                hough_space(rho_index, theta_index) = ...
                     hough_space(rho_index, theta_index) + 1;
            end
        end
    end
end

%% show hough transform
subplot(1,2,1);
imagesc(theta, rho, hough_space);
title('Hough Transform');
xlabel('Theta (radians)');
ylabel('Rho (pixels)');
colormap('gray');

%% detect peaks in hough transform
r = [];
c = [];
[max_in_col, row_number] = max(hough_space);
[rows, cols] = size(image);
difference = 25;
thresh = max(max(hough_space)) - difference;
for i = 1:size(max_in_col, 2)
   if max_in_col(i) > thresh
       c(end + 1) = i;
       r(end + 1) = row_number(i);
   end
end

%% plot all the detected peaks on hough transform image
hold on;
plot(theta(c), rho(r),'rx');
hold off;

%% plot the detected line superimposed on the original image
subplot(1,2,2)
imagesc(image);
colormap(gray);
hold on;

for i = 1:size(c,2)
    th = theta(c(i));
    rh = rho(r(i));
    m = -(cos(th)/sin(th));
    b = rh/sin(th);
    x = 1:cols;
    plot(x, m*x+b);
    hold on;
end

  Linked:
  How to do De-Houghing of a Hough transform'ed Image?

Firstly, Matlab has a Built in Hough Transform: no need to reinvent the wheel.

[H,T,R] = hough(BW,'RhoResolution',0.5,'Theta',-90:0.5:89.5);

Although your image doesn't necessarily require edge detection you could improve the processing time and effectiveness of the algorithm my using it. Your triangle has fat regions of white and black. Ideally, the triangle would be 1px thick marking the edges of the triangle. Use Canny Edge Detection

BW = edge(Image,'canny');

the result is \theta is in the range $-90 < \theta < 90$ degrees while your plot is $0 < \theta < 180$ (or $0 < \theta < \pi$) so subtract $90$ degrees ($\pi/2$) .

There is the potential for you to choose the wrong peak because there are neighboring peaks that may be larger in the accumulator matrix. While there are many algorithms here is one that I've seen used in Hough Transforms in the past:

1) Define a region shape (typically its square)
2) Define an accumulator threshold
3) Select one pixel in the accumulator matrix
4) If it is a peak (i.e., larger than neighboring values above a threshold)
       a) look at the points within the region shape.
       b) If there is a larger value
              this is not a peak
          Else
              this is a peak
 5) Move to next pixel in accumulator matrix.

Look into HoughLines for displaying the hough transform lines, results:

http://www.mathworks.com/help/toolbox/images/ref/houghlines.html

Effects of Using Canny Edge Detector

Edge Detection can potentially turn each side of the triangle into two lines.

The goal of canny edge detection is to produce maximally thin/narrow edges by using nonmaximal supression

Canny Edge Detection in a Nutshell (Source: Digital Image Processing, Gonazalez)

1) Smooth input Image using a Gaussian Filter
2) Compute the Gradient magnitude and angle (Sobel, Perwitt or robert cross filters)
3) Apply Nonmaxima suppression (this is where the thinning happens)
   a) Figure out which direction the edge is
   b) If the edge's magnitude is smaller than one of its two neighbors in the direction of the edge
          set the edge point to zero
      Else
          leave it alone
4) Use double thresholding and connectivity analysis to detect and link edges


Designing complex magnitude compensating FIR filter I am trying to find a Python replacement for the Matlab function cfirpm (complex FIR Parks-McClellan). I use cfirpm in the design of a complex FIR magnitude and group delay compensation filter for a crystal IF filter. Since the frequency response of the crystal IF filter is not symmetric, the compensating filter isn't either. Therefore, I cannot simply use a lowpass filter design and multiply it by e^(jwft). I have searched for hints on how to do this on this site and many others. I have looked at 23147/design-filter-with-arbitrary-magnitude-and-phase-response-in-matlab where the answer was to find the inverse FFT of the two-sided complex frequency response, etc. but what I have from a frequency sweep using an RF signal generator and processing the I and Q samples is a vector representing the desired magnitude at evenly spaced frequencies and I'm not sure how to represent that as a two-sided complex frequency response. Any ideas on how this type of filter can be designed in Python would be greatly appreciated.

Here is an example I already had in Python demonstrating the Wiener-Hopf equations detailed here Compensating Loudspeaker frequency response in an audio signal to solve for the channel equalizer. This application by the OP is ideal for this given the static channel condition as long as deep spectral nulls don't exist, over a iterative LMS or RLS equalizer solution such as those detailed here (which are great for dynamic channels that don't exhibit deep frequency nulls in band): Recursive Least Square Adaptive Linear Equalizer
For channels with deep spectral nulls in band (frequency selective fading), a decision-feedback equalizer would be the better approach due to degradations from noise enhancement with the linear equalizer approaches.
Eye diagram of undistorted waveform with no channel distortion, as it would be received after the 2nd RRC filter in the receiver:

Tx Spectrum of undistorted waveform prior to channel distortion:

Example complex channel distortion that we will equalize:
cir = [.019*np.exp(1j*.2), .021*np.exp(-1j*.24), .08*np.exp(1j*2.4),  .032*np.exp(-1j*.14), -.032*np.exp(-1j*.14)]

This is the same eye diagram and Rx spectrum above after passing through our channel (yuck, what a mess! How do we possibly demodulate this????). Notice the assymetrical passband, a tell-tale sign that the compensating filter MUST be complex.

Compute Equalizer Solution Using Tx and Rx Waveforms
(Note this is Python 3.8, but nothing jumps out at me that wouldn't be compatible with earlier versions)
import scipy.linalg as linalg
import numpy as np

def convmtx(h,n):
    # creates the convolution (Toeplitz) matrix, which transforms convolution to matrix multiplication
    # h is input array,
    # n is length of array to convolve
    return linalg.toeplitz(np.hstack([h, np.zeros(n-1)]), np.hstack([h[0], np.zeros(n-1)]))

# Compute equalizer given known tx and rx waveforms

omit = 100    #initial samples to exclude
shift = 0  # number of samples to shift dominant equalizer tap to the right

ntaps = 20
depth = ntaps * 30
delay = ntaps//2
A = convmtx(rx[omit+shift: omit+shift+depth], ntaps)
R = np.dot(np.conj(A).T, A)
X = np.concatenate([np.zeros(delay), tx[omit: omit+depth], np.zeros(np.int(np.ceil(ntaps/2)-1))])
ro = np.dot(np.conj(A).T, X)

equal_coeff = np.dot(linalg.inv(R),ro)

Below is the resulting frequency response for the equalizer and the equalized waveform. The details for the parameters used such as equalizer number of taps and depth are further detailed in the earlier links given. Note how the equalizer can only determine the compensation where signal energy exists (otherwise it has no way to measure the channel), therefore it is important to use a test (sounding)  waveform that fully covers the bandwidth of the channel of interest. Pseudo-random sources or any white noise source is ideal given the broad spectral coverage.


correlation matrix for complex white noise does not equal to identity matrix? I generate 8 time series of gaussian white noise and want to compute the correlation matrix at a frequency bin between them. The following is the matlab code I use
clear all;
close all;

fs = 4000;
N = 100*fs;
R = 0.5;% overlap percent
Nseg = fs;

M = 8;

for m = 1:M
   x(:,m) = randn(N,1) ;
end

segNum = floor((N-R*Nseg)/((1-R)*Nseg));

for n = 1:segNum
   y = x([1:Nseg]+(n-1)*Nseg*(1-R),:) ;

   Y = fft(y);

   RR(n,:,:) = transpose(Y(300,:))*conj(Y(300,:));

end

R_hat = 0;
for n = 1:100
   R_hat = R_hat + squeeze(RR(n,:,:));
end

R_hat = R_hat/100;

However, the estimated correlation matrix $R_{hat}$ does not equal to identity matrix. Can you help point out what the problem is with my code?

Why don't you simply compute $X^{H}  X / N$ ?
clear all;
close all;

fs = 4000;
N = 100*fs;
R = 0.5;% overlap percent
Nseg = fs;

M = 8;

I noticed your data was not complex, since you metioned complex
in the question here I made it complex.
X = 0.5*randn(N,M) + 0.5j *randn(N,M);

This will be close to the identity matrix there will be small complex perturbations
the expected perturbations are inversely proportional to $\sqrt{N}$. The $X^H$ corresponds to X', it is the conjugate transpose of $X$.
R_hat = X' * X / N

Edit
Here I compute RR such that RR(:,:,k) is the accumulated covariance matrix of the frequency bin $k$ for the different inputs.
clear all;
close all;

fs = 4000;
N = 100*fs;
R = 0.5;% overlap percent
Nseg = 256;

M = 8;

for m = 1:M
   x(:,m) = randn(N,1) ;
end

segNum = floor((N-R*Nseg)/((1-R)*Nseg));

RR = zeros(M, M, Nseg);

for n = 1:segNum
   y = x([1:Nseg]+(n-1)*Nseg*(1-R),:) ;
   Y = fft(y);
   for k = 1:Nseg
     %% the correlation for for one bin one slice.
     RR(:,:,k) = RR(:,:,k) + Y(k, :)' * Y(k, :);
   end
end
%% Normalize
RR = RR * (1/segNum/Nseg);

%% Acculate over all bins
R_hat = 0;
for n = 1:Nseg
   R_hat = R_hat + RR(:,:,k);
end

R_hat = R_hat/Nseg


Constructing a Gaussian kernel in the frequency domain I'm currently learning about Fourier transform, but find the differences between spatial domain and frequency domain a bit confusing at times.
Let's say I would like to perform convolution of an image with a Gaussian kernel. As far as i understand the Fourier transform of a Gaussian is also a Gaussian i.e the Fourier transform of
$$
g(x,y;\sigma) = \frac{1}{2\pi\sigma^2} \cdot e^{- \frac{x^2+y^2}{2\sigma^2}}
$$
is
$$
G(x,y;\sigma) = e^{- 2\pi\sigma^2(x^2+y^2)}
$$
However if I directly construct a kernel using $G(x,y;\sigma)$ the output kernel looks completely different from first computing $g(x,y;\sigma)$ and then doing the FFT. What is the correct approach here? I hope the question is clear, but it is essentially how I can construct a gaussian kernel directly in frequency space so i do not have to FFT the kernel.

You can check numerically that fft(g) = G, that will be easier to analyze than the filtered result. At least you can easily plot it and identify if the coefficients are shifted, flipped, or scaled. One thing easy to do incorrectly is the determination of the frequency corresponding to each coefficient. I would strongly recommend that you first aproach the 1D version of your problem, here you have an example of 1D filtering I wrote for a related question.


MMSE Equalizer in communication system This question is regarding applying minimum mean square operation at the receiver side of communication system.

I have read in paper that

Given an output signal from communication system $$R(k)= H(k)S(k) + N(k)$$

an MMSE filter $$G(k)= \frac{H^*(k)}{|H(k)|^2+\frac{\sigma_n^2}{\sigma_s^2}}$$

My incomplete solution

I think we need to find G(K) such that

$$\mathbb{E}[G(k)H(k)S(k)+G(k)N(k)- S(k)]^2$$

Does anyone know how should I continue, thanks alot.

For simplicity of notation I'll leave out the index $k$. You've already noted that the error is

$$E=R\cdot G-S=S(HG-1)+NG\tag{1}$$

The MSE can be computed by integrating the power spectrum of the error. Since the power spectrum is real-valued and non-negative, its integral (i.e., the MSE) can be minimized by making the power spectrum as small as possible. Assuming white data and white noise with powers $\sigma_s^2$ and $\sigma_n^2$, respectively, which are independent of each other, we can compute the error power spectrum from (1) as

$$S_E=\sigma_s^2|HG-1|^2+\sigma_n^2|G|^2\tag{2}$$

The expression for the equalizer response $G$ that minimizes (2) can be found in an elegant way by completing the square in (2):

$$S_E=(\sigma_s^2|H|^2+\sigma_n^2)\left|G-\frac{H^*}{|H|^2+\frac{\sigma_n^2}{\sigma_s^2}}\right|^2+\frac{\sigma_n^2}{|H|^2+\frac{\sigma_n^2}{\sigma_s^2}}\tag{3}$$

The expression (3) can be minimized by choosing $G$ such that the left-hand term vanishes:

$$G_{opt}=\frac{H^*}{|H|^2+\frac{\sigma_n^2}{\sigma_s^2}}\tag{4}$$

which is the optimal MMSE solution.

It is not immediately obvious that (2) and (3) are identical, but proving their equivalence is simple (yet tedious) by multiplying both expressions out and comparing them.


Find the valley in the given in a sequence of numbers I am trying to find the valley in a sequence of numbers of around 2000 length.

Always prominently there will be two peaks. But the value and location might change. I need to find the valley, defined as the lowest point between the two peak values.

I tried using the following code:

% sums = list of positive integers generated from an analysis
inverted_sums = double(max(sums) - sums); % inverting the sequnce
[~,locs] = findpeaks(inverted_sums) % finding the peaks

But the above code returns many valleys with even a small difference, including the correct one. How do I solve it?

The following code will find the lowest point between two prominent peaks.

%Smoothing the data
movAvgFiltOrder = 8;
b = filter(ones(1,movAvgFiltOrder)/movAvgFiltOrder,1,sums);

% finding peaks from the smoothed curve
[vals,pos] = findpeaks(b);
idx = find(vals>mean([min(vals),max(vals)])); % removing peaks which are below the average,
%this is helpfull if some lower amplitude peaks are caught at sides of the prominent peaks.

peak_pos1 = pos(idx(1)); % farthest prominent peak on one side of valley
peak_pos2 = pos(idx(end)); % farthest prominent peak on other side of valley

% correcting the peak positions for filter delay (Coarse tuning)
peak_pos1 = peak_pos1 -  movAvgFiltOrder/2;
peak_pos2 = peak_pos2 -  movAvgFiltOrder/2;

% Looking for any other peaks nearby (Fine Tuning)
tuning_range = min(movAvgFiltOrder,floor(abs((peak_pos2-peak_pos1-1)/2)));
p1 = peak_pos1-tuning_range:peak_pos1+tuning_range;
p2 = peak_pos2-tuning_range:peak_pos2+tuning_range;
[~,r1] = max(sums(p1));
[~,r2] = max(sums(p2));
peak_pos1 = p1(r1);
peak_pos2 = p2(r2);

% Identifying valley
valley = sums(peak_pos1:peak_pos2);

% finding minimum point in valley
[min_val,pos_min_val] = min(valley);
valley_amplitude = min_val; % amplitude of minimum point
valley_position  = peak_pos1+pos_min_val-1; % finding the valley point in the entire array

This valley_position will give the poisition of the lowest point between major peaks.


How many coefficient sections should I select when designing filters with MATLAB FDATOOL? I attach a picture about the results of a filter design in Matlab.

It says that it has 12 sections. One section contains the coefficients of the transfer function. If I'm implementing this filter in a microcontroller, I think that 12 sections would take much time as for to do this inside an interrupt loop.

What happens if I only choose 3 sections? Would this affect badly the filter performance?

If you had infinite-precision math, then no it would not change anything to split your order-24 IIR in 3 cascaded order-8 sections. You could even implement it directly as a single difference equation, i.e. no cascade.

However, in the real world, we don't have infinite-precision math, we typically use 64-bit doubles, 32-bit singles or even integers (fixed-point). In those cases, implementing an order-24 IIR as a single difference equation would lead to quantization issues. Your filter would likely be unstable. Even if it were stable, it would probably not have the expected frequency response.

Therefore what we usually do is split the order-2M IIR filter into M cascaded order-2  IIR filter called byquads. The quantization issue is much easier to deal with when using byquads, at least when using doubles or singles. With Fixed-point arithmetic, even when using byquads, you need to double-check everything.

As a bonus, in DSP libraries, you usually have optimized functions for byquads IIR.

Care to explain why you need an order-24 IIR? I fail to see why you would need something more than an order-8.


Calculating SNR from Frequency-time domain I tried to analyse the time-frequency of specific appliances. As it shows in the following frequency-time graph
I have High frequency noise from different appliances here. 4096 is the number of frequency divisions. So in the data I have, across a spectrum of 0-1Mhz, there are 4096 divisions. In other words, each "bin", holds the amount of energy at that frequency. 1Mhz/4096 = 244.14Hz, thus, each bin from 0-4095 represents energy averaged over roughly 244.14Hz. (244,14 vectors per second)

I tried to analyse each bin before the time before 1.34334165, and plotted the 4096 point FFT vector as frequency vs. power. I am not sure why it is 500 Khz, but I think the High Frequency signal is sampled at 1 MHz and FFTs are computed to get 2048 frequency data points from DC to 500 kHz

So the question is, I need to calculate the SNR of each bin represented in the graph. I was following the thread here on how to calculate SNR from frequency domain signal, but I was wondering if this answer is correct. I couldn't understand the 50 leakage and, as I was wondering too, is the inputsignal the  sum of the magnitudes of all the 4096 bins?

This is probably not going a particularly satisfactory or useful answer for you. SNR is a ubiquitous concept in Signal Processing and one would expect a correct straightforward calculation on data would be available and it basically is:
$$
SNR=\frac{\text{Signal Power or Energy}}{\text{Noise Power or Energy}}
$$
so you just look at your spectra and identify the signal and noise and make the calculations. This is usually straightforward to do theoretically. In some applications these identifications are easy and in others less so and looking at your data, not much looks obvious. A FFT of the data didn't make what is signal or noise obvious to me and if were obvious to you, you probably would not be asking so
$$
\text{Correct} \;SNR \ne \frac{\text{Hard to tell Signal Power or Energy}}{\text{Hard to tell Noise Power or Energy}}.
$$

Instead,
$$
\text{Plausible} \;SNR = \frac{\text{Plausible Signal Power or Energy}}{\text{Plausible Noise Power or Energy}}.
$$
is likely to be more successful. I've actually had to present these issues to management and at the end been asked if this was "correct". Saying yes has usually been the "correct" answer.

Getting back to you data, the frequency time graph has some quiet spots that are probably your best candidates for a plausible noise estimate.  This assumes, that the noise is stationary, but looking at everything else in the time frequency plot, this seems plausible. Extract a Noise Power estimate from the quiet segments and average as many as possible. Also I'm guessing that the signals are of finite duration.
$$
\text{Plausible} \;(S+N)R = \frac{\int_T (\text{Plausible Signal + Noise Power}) \; dt}{\text{Plausible Noise}\times T}.
$$
where $T$ is the signal duration.

Some details, We have Parseval's Identity,
$$
\sum_{n=0}^{N-1} x[n] = \frac{1}{N}\sum_{k=0}^{N-1} \mid X[k] \mid ^2$$
and the variance of a sum of $N$ independent random variables
$$
\sigma^2_{\text{Total}}= \sum _{k=0}^{N} \sigma_k^2.
$$
The random variables don't have to be Gaussian, they just need to have a variance.

For long DFT's noise bins are approximately independent, summing the noise power present in each bin is a reasonably same assumption if the noise is approximately flat across the spectrum.

If your signal is known, the energy can be calculated using Parseval's Identity.  If your signal has random characteristics, The signal bins are probably not independent and summing bin power may not be the best plausible method.  If you have enough data, you can estimate the covariance
$$
Pwr=\mathbf{x}_{sig}^H \mathbf{R}\mathbf{x}_{sig}
$$
which is a power quantity.

For complicated signal and noise scenarios the most reasonable method is to derive a model, not necessarily a full model, and simulate your SNR estimation process.

The answer to 50 is a heuristic to account foe signal leakage outside the nominal signal limits.

Two issues with your question

The time frequency plot indicates frequencies to 1MHz but your spectrum is 500KHz.  You really should know why.

Using Unix time as a time index is terrible.  You make the reader work.

Most of this didn't really give you the answer you were looking for.  There is no straightforward fully general way to calculate SNR without knowing what is signal and what is noise to a high degree of specificity.


Is ergodic in mean a property defined only for WSS stochastic processes? I understand the definition of a random process $X(t)$ being ergodic in mean (first-order ergodic) is that the expectation of the sample mean $<u_X>_T=\frac{1}{T} \int_{\frac{-T}{2}}^{\frac{T}{2}} X_w(t)dt $ should converge to the ensemble mean as $T \rightarrow \infty$. Can we define ergodic in mean property for processes which do not have constant mean.

I know that the property of being ergodic in ACF, second-order ergodic, requires WSS assumption.

My question is motivated by the fact that the requirement of WSS is not given in the wiki article describing ergodic process but used in the example illustrating the same.

Time average is defined over a time interval, T which will always result in a constant value. Keeping this constraint in mind, ensemble mean is always forced to be constant for satisfying the condition of ergodicity.


How to implement a gradient based Hough transform I am trying to use the Hough transform for edge detection, and would like to use gradient images as the basis.

What I have done so far, given the image I of size [M,N] and its partial derivatives gx, gy, is to calculate the gradient angle in each pixel as thetas = atan(gy(x,y) ./ gx. Similarly I calculate the gradient magnitude as magnitudes = sqrt(gx.^2+gy.^2).

To build the Hough transform, I use the following MATLAB code:

max_rho = ceil(sqrt(M^2 + N^2));
hough = zeros(2*max_rho, 101);
for x=1:M
    for y=1:N
        theta = thetas(x,y);
        rho = x*cos(theta) + y*sin(theta);

        rho_idx = round(rho)+max_rho;
        theta_idx = floor((theta + pi/2) / pi * 100) + 1;
        hough(rho_idx, theta_idx) = hough(rho_idx, theta_idx) + magnitudes(x,y);
    end
end

The resulting Hough transform looks plausible (see http://i.stack.imgur.com/hC9mP.png), but when I try to use its maxima as edge parameters in the original image, the results look more or less random. Did I do something wrong in constructing the Hough transform?

UPDATE: I had a stupid mistake in my code: rho was calculated as x*cos(theta)+y*cos(theta) instead of x*cos(theta)+y*sin(theta). That is, I was using two cosines instead of a cosine and a sine. I have edited the code above and the new resulting image is below. This did not give much better edges though.

@endolith: To plot an edge, given a maximal value in the hough-matrix at rho_idx, theta_idx, I translate the indices to rho,theta values:

theta = (theta_idx -1) / 100 * pi - pi / 2;
rho = rho_idx - max_rho;

Finally I plot the edge as y= (rho - x*cos(theta)) / sin(theta).

I am a bit confused by your question.  Hough transform is used to detect lines, not edges.

If all you want is an edge map, you should simply threshold the gradient magnitude, or use something more fancy like the Canny edge detector.

If you want to detect straight lines, you would be better off starting with an edge map, and then using the hough function if the Image Processing toolbox, if you have access to it. The problem with doing a Hough transform on the gradient is that edge pixels forming a straight line might have opposite gradient orientations.  For example, consider a checkerboard pattern.  A edge between two rows of squares flips orientation depending on whether you have a black square above and a white square below, or the other way around.

As to your implementation, I think the problem is that the bins in your Hough matrix are too small.  Essentially the bin size in the rho dimension is 1, and the bin size in the theta dimension is less than 2 degrees.  That means that the gradient orientations must alight very precisely to form a line, which rarely happens in practice.  If you caculate rho_idx and theta_idx so that the bins are larger, that will make your line detector more tolerant to errors, and you might get better lines.


How to count the number of threads in a bolt or nut using image processing? How can I count the number of threads in a Nut with image processing?

I know that I need to find the place of threads and count them. Ok, But how can I find the place of threads?

This is a sample image:

Update
As I know I need to do these steps:

Making image ready. For example az @A_A sayed, with Thresholding
Remove noise and just keep threads in the image
Count threads. maybe using signal processing

Level 1 is solved. What should I do for 2 and 3?

For example this is what I did for a bolt. I got a sample of a vertical line in the middle of the nut then plotted the data. But I don't know how to remove additional parts and count up and downs.

Thanks

The most important thing in this problem is to be able to control the scene and more importantly:

The lighting
The camera position & angle

The Lighting:
Simply put, the lighting is what creates the essential contrast to delineate the thread. If you light the thread "head on" (light from the viewpoint of the camera), there will be no contrast or "ups and downs" to count. The better position for the light is at the top and slightly forward of the face of the thread so that a) All thread peaks are lit (both those that are on the top, but also those that are below them) and b) The top thread casts a bit of shadow to the one below it. Try to think of the thread as a difraction grating, only here, you need the shadows.

The Camera: The camera distance from the thread and angle to the thread are very important, especially if you wish at some point to derive the pitch of the thread. To derive the pitch of the thread, you would have to count lines per unit of distance.

The Camera Distance: The only way to relate the pixel size to a physical unit of length in a monocular image, is to fix the distance between the subject and the lens, similarly to how they do it in archaeology photography. There, we assume that the ruler's plain coincides with the focus plane of the lens, or is too near to be creating a negligible error.
The Camera Angle: The camera angle is incredibly important in this problem because it is responsible for an optical distortion that would affect your measurement. In the nut example, it is impossible to image the thread vertically ("head on"), because you are trying to image the inside of a very small and thick cylinder. Therefore, you have to look at it from an angle. In the bolt example, it is impossible to image the thread vertically ("head on"), because there will be no shadows to create contrast. A pen, imaged at an angle will appear shorter than it really is. Similarly, the thread imaged at an angle, will not give you equidistant peaks. Specifically, as the thread "screws away" from the camera, its pitch will be increasing because the lines that denote the threads would appear to get closer and closer and closer. Therefore, the angle that the camera is set against the screw determines a) The error of measurement and b) the limits of measurement (how fine / coarse a thread you can measure).

AFTER determining / fixing these parameters which are very important to the problem you are dealing with, you can start thinking about extracting the thread. If you get the imaging conditions right, the thread extraction via image processing will be much more easier and straightforward than trying to do it over any random image.

Trivial way of extracting the thread:

Assuming that the thread is aligned and centred within the camera's Field Of View (FOV), you are using something like Octave and the image has been loaded in I with something like I = imread('someNutImage.jpg');:

Threshold the image

Because the thread is (supposedly) aligned with the FOV, take the middle column of values. That would look something like:

Notice here how the pitch changes visibly because of the camera angle.

Now, "count" the gaps. This is a trivial "pulse counting" excercise which I have outlined previously. Prior to counting the pulse widths, you need an absolutely saturated signal. Therefore Q = double(I(:,320)>128).*255; (Where 320 is some column image, I am taking the middle here as explained above). Now, differentiate the signal to get the "pulse fronts" and then differentiate that again, to get the "inter-pulse-front" distances. That is:

pulse_fronts = find(abs(diff(Q)));
pulse_lengths = diff(pulse_fronts);

For the signal I am using here, I get:

1, 24, 64, 11, 8, 13, 6, 14, 6, 13, 6, 13, 6, 13, 5, 12, 6, 11, 5, 12, 5, 10, 6

That is, the signal starts with a spike pulse of length 1, followed by 24 zeroes followed by 11 ones, followed by 8 zeros, followed by 13 ones and so on.

Of course, you are not going to "count the gaps" manually. At this point, you can do a histogram and obtain the two most frequent values. For example, with a hist(pulseLength,max(pulseLength)), you get something like:

The "low" most frequent value (6) are the zeroes and the "high" most frequent value (13) are the ones. Once you have these two values (from the histogram), you know the full length of your thread in pixels. To convert it to mm, you need to determine the pixel size.

Less trivial way of extracting the thread:

The above method allows you to filter (if required) the pulse lengths to get rid of various artifacts. Another way to extract the thread would be to use the Discrete Fourier Transform of Q to assess the frequency of the pulse train.

For example: semilogy(abs(fft(Q))), gives you something like:

This image is 480 pixels tall. This means that the highest frequency (in terms of alternating black/white lines per unit of length) is 240. The first peak of the DFT occurs at 26 lines. From the method above we had ~13+~6=~19. The DFT will give you an average pitch over the whole acquired signal. Again, to give these numbers a physical meaning, you would need to know the pixel size which is only known by fixing the distance from the camera.

Even less trivial way of extracting the thread:

Once you threshold the image, either apply thining or a simple high-pass convolution mask to preserve the edges and then a Hough Transform. This will generate one "dot" for every slanted "line" that is created by the thread. What you would be looking for then would be a series of dots, one next to the other signifying the thread. By measuring the distance between these dots, you can extract your thread information. This is probably the most accurate method because the Hough Transform is very sensitive and accurate at recovering lines from very noisy images like the thresholded one we are talking about here. It is also the most computationally expensive option.

Each one of these methods has pros and cons for estimating the thread of a nut or bolt.

Finally, the most important part in such applications is to fix the imaging geometry so that you know what you measure and what sort of limits this measurement has.

Hope this helps.


What exactly is "sparse representation"? I saw a recommended topic for the final project in my university called "dsp and dip applications using sparse representation techniques (MATLAB, C, C++)".

I consider  taking this topic as my final project, but there are not so many things about it on the internet to understand what this is exactly.

So, I would like to hear some clear answers about it and some opinions about how productive would that be for me considering that I want to learn more and more about dsp and dip.

Take a sine-like signal $s$. In the appropriate Fourier $\mathcal{F}$ domain, it is represented by two "peaks", the other coefficients being zero. Fourier is a sparse representation for sines or close-to-sine signals. Conversely, a zero signal, except for a few values, is sparse in its original domain.

In narrow sense, a sparse representation of data is a representation in which few parameters or coefficients are not zero, and many are (strictly) zero. This can be measured by the $\ell_0$ count index, which yields the number of non-zero components. Here, $\ell_0(\mathcal{F(s)})=2$.

However, in practice, it is very rare to have few exactly non-zero components, and a lot being exactly zero. This happens especially as floating point computations are rarely exactly zeros, albethey very small. You can see for instance the simulations for Why is $\exp(\ln(x))-x \neq 0$ in floating point arithmetic?

So, broader senses are accepted in practice, amounting to few high components, many small components.
I will thus extend the discussion on terminology on either "sparse" and "representation".
On sparsity, or parsimony, one should mention the economy principle, law of parsimony or  Ockham's razor:

  An explanation of a thing or event is made with the fewest possible
  assumptions

Initially, this pre-scientific principle has being used in theological debates.
In science, it is  used as a heuristic guide in the selection or development  of theoretical models. One early example is the geocentric model. Starting from the belief that Earth was the center of the Universe, further observations revealed complicated star motion patterns, building a model with complex (epicycloidal) curves, requiring many parameters to explain. Until some understood that the hypothesis of a centered earth was not appropriate, and that other models could provide simpler explanations, and formulae, with the same, or a better, prediction or description power. Kepler or Newton laws were way more efficient.

Another example is polynomial representation: suppose that you have fifteen couples of $x_i,y_i$ values that you want to use for prediction. Of course a Lagrange polynomial of degree 14 will fit then perfectly. Should it be used for prediction? Probably not, as data are imprecise, and the polynomial extrapolation will be unstable; maybe, a one or two degrees fit, according to some cost function, might be more appropriate. The idea is that two or three parameters (slope/intercept/second degree term) could be more efficient for a given task.

We now have the major ingredients. Take a set of data $d_i$,   with meaningful content. The data generally has many dimensions, $i\in \mathcal{I}$ (think about the number of pixels in images). There is some expectation that there exists a combination of a small subset of features (sparsity) $c_k$ that can faithfully describe  the dataset, or predict, infer from this dataset. In other words, with $k\in \mathcal{K}$, with cardinals such that $|\mathcal{K}| \ll |\mathcal{I}|$, there is a transformation $T$ such that the quantities $d_i - T(c_k)$ are small enough.

In the case of sines, $|\mathcal{K}|=2$ is sufficient in a Fourier basis. In practice, one tries to  have sequences of $c_k$ such that most of them are very small (if not zero), and few of them recovers data correctly. There are measures to quantify the sparseness of a sequence, and to quantify how close the sparse representation is to the actual goal.

Transformations  to represent data are standard tools from harmonic analysis (Fourier, Wavelets, filter banks, Radon, bases or redundant frames), statistical analysis (least-square regression, PCA, PLS, lasso, ridge regression), clustering and segmentation tools, etc. The distance is often  some norm or quasi-norm. Even neural networks strive to provide learned sparse representations from a training set.

Sparsity indices are 0-homogeneous (Gini index, norm ratios) or akin to norms or quasi-norms for minimization purposes. Some use a milder description, i.e. the compressibility, reordering coefficients $T(c_k)$ with decreasing magnitudes $|(c_{\sigma(k)}|$ and studying the exponent $\alpha\le 0$ such that:

$$\|d_i - T(c_{\sigma(k)})\| \leq K|k|^{\alpha}$$

or "how fast   compressible representations can approximate a data or a dataset".

I consider the idea of sparse representations to be a key driver in machine learning, classification, compresison, and you can find a plenty of ressources to reuse.