Find horizontally and vertically aligned sets of points I have the image below.  This image contains sets of points that are aligned horizontally and vertically.

How can i find these sets of points? based 1) on their alignment and  2) on the ** **inter-line and Inter-word distances.

Do I have to use these 2 arguments, alignment and distance? or is one of them sufficient?

Is it possible to use the Hough transform? or another algorithm?

Using alignment, it's mostly a matter of heuristic (sweeping horizontally and vertically). You can also work out something similar to Hough transform here.

Using inter and intra-dot distances is also known as linkage. What you want to do can be obtained through hierarchical clustering using a single linkage distance computation. It works as follows:

you start with all your dots as (unconnected) leaves in a tree;
you test all the leaves pairs, compute their spatial distance, and connect the two leaves that yield the smallest distance;
replace these 2 leaves by only 1 leave that represents the cluster;
go back to 2 and repeat until you have only 1 leave left.

The last leave is the root of a tree that represents the hierarchical clustering result, which is also called a dendrogram.
By using single linkage, you connect each dot to the closest dot, then recursively all closest pairs of dots, etc. You can finally discard the hierarchical tree part that corresponds to the smallest distance between two different lines.


Parallel Cascaded FIR Filters I have a relatively simple question. I know if you connect multiple LTI systems in parallel, overall system's frequency response is the sum of each LTI system's frequency response. Therefore, what would be the frequency response of the overall system if two FIR Low-Pass filters are connected in parallel? Let's say one has a cut-off frequency of 50 Hz and the other one has 100 Hz. I'm guessing it would have a frequency response of 1 between 50-100 Hz and 2 below 50 Hz.

Similarly, how would you construct a band-pass filter from a pair of low-pass filters? I am totally stumped with this one.

As you've mentioned, you can simply add the impulse responses, and, consequently, the frequency responses of LTI systems connected in parallel. In general, frequency responses are complex-valued, so you can't just add their magnitudes. However, if their phase responses are identical, you can add their magnitudes. A special case are two linear-phase FIR filters of the same length, which have the same (linear) phase response, so their magnitudes can be added.

If both filters approximate a magnitude of $1$ in their passbands, and a magnitude of zero in their stopbands, then the resulting filter will approximate a magnitude of $2$ where their passbands overlap, a magnitude of $1$ where the passband of one filter overlaps with the stopband of the other, and a magnitude of zero where both stopbands overlap. Note that FIR filters cannot have a constant magnitude (unless they have only $1$ filter tap), so all constant values are just approximated.

Now that you've figured out the resulting response when the outputs of the two filters are added, maybe you can come up with a creative way to connect the outputs of two linear phase FIR low pass filters (of the same length) with different cut-off frequencies to implement a bandpass filter.


Correlations of matrix I have a matrix of 10 column, each column has measurements at different times. I have a second matrix also of 10 columns, each consisting of measurements at different times.
I want to carry out cross correlations of every combination within this matrix. Example, xcorr row 1 matrix 1 with row 1 of matrix 2, then row 1 of matrix 1 of row 2 of matrix 2 and so on. I would like correlations and time lag for each and every combination of columns.
I currently do this, but it means I have to do each row separately which is very time consuming.

[xc,lags] = xcorr(matrix1(:,1),matrix(:,2));
[~,I] = max(abs(xc));figure;
plot(lags,xc);
legend(sprintf('Time delay %d',lags(I)));
title('Column 1 matrix 1 xcorr with Column 2 matrix 2')

matrix1 = rand(10, 100);
matrix2 = rand(10, 100);
for ind1=1:size(matrix1, 2)
    for ind2=1:size(matrix2, 2)
        [xc{ind1, ind2},lags{ind1, ind2}] = xcorr(matrix1(:,1),matrix2(:,2));
        % Here you can add any other operation you want
    end
end


How to perform model fitting for system identification I am having a really hard time in understanding how to formulate a model say linear AR model to represent a communication channel or maybe any motion. I have the experimental data representing the kinematics of robot movement. In the case of a communication channel, it will have a medium, transmitter and receiver. So, will the medium be represented by AR model or do we represent the transmitter and receiver by these models? How to decide the order of the system to be chosen initially. I am not from the area of signal processing but need to use this for a part of my research work. I have gone through literature review and the background of Kalman filter, Least square method, recursive least square for estimation purpose. But I cannot find a good example which shows how the parameters are estimated from the time series or motion data or speech modeling. Can somebody explain or point out the starting point and how to go about it? I know this is a very broad topic but an explaination with a small example will be really an eye opener for a beginner.

In any signal processing problem, there are usually two components: the signal model and the channel model.

The signal model is the mathematical description of how your ideal signal, call it $s(t)$ is generated.

The channel model is the mathematical description of how your channel corrupts or alters the signal.

One aim in telecommunications is to look at the received signal $r(t)$ --- $s(t)$ corrupted by the channel --- and try to undo the effects of the channel.

For example, suppose
$$
s(t) = A \sin(\omega_0 t + \phi)
$$
and suppose the channel is:
$$
r(t) = h(t)*s(t) + n(t)
$$
where $h(t)$ is the impulse response of the channel, and $n(t)$ is some additive noise.

In this context, your "AR"  filter is the $h(t)$.

If the context is in robotics or kinematics, then things get a little more complicated.

Did you open the "Algorithm" tab on the link in your comment? It's pretty straighforward:

  The polyfit MATLAB file forms the Vandermonde matrix, $\bf V$, whose elements are powers of $x$:
  $v_{i,j} = x_i^{n - j}$
  It then uses the backslash operator, $\backslash $, to solve the least squares problem ${\bf V} p \approx y$.


When can a windowed cosine be considered band-limited signal? I have an exam in the following days and i have no clue of what is the response to this question in the exam:

  given this signal:

  $$g(t) = H\left(t+\frac{T}{3}\right) f(t) - H\left(t-\frac{2T}{3}\right)f(t)$$ $H(t)=$ Heaviside step function

  $f(t) = \cos(\omega_1 t)$

  Define a value of $T$ for which the signal $g(t)$ can be considered
  band limited.

This question doesn't make much sense to me. This is a windowed cosine function and it's bandwith depends on the value of $T$, it's always limited!
Unless $T$ approches 0 and it can be seen as an impulse. Am I missing something?

T=0 or T=infinity.

Any non-zero signal with finite length support in one domain has infinite support in the other domain.

cos(wt) is only theoretically band-limited if infinite in length.


Use of I/Q representation for unmodulated Rx Signal it may be a stupid question but while learning about DSP and Modulation I was wondering, if the I/Q representation of an RF Signal is of any use, when the received signal does not contain any modulation.
So I understand the concepts of Inphase and Quadrature components, how they are modulated in the sender, combined for transmission and demodulated in the receiver. But if i where to receive some other signal, lets say for example RF-emissions from some random or physical source i want to observe. Does the separation into I- and Q components yield any value then? Since it has never been composed to hold information by the source. Or would there still be a value in having I- and Q components?
I couldnt find any explicit statement on that yet.
Cheers
Lucas

If you are "observing" the source, this implies there is some sort of information you are looking to get out of it, whether it be the total background noise, interference levels etc.
Do you find "value and use" in the RF signal's magnitude versus time? What about the RF signal's phase versus time? The IQ representation gives us both of these and it simply is another way to represent complex numbers:
$$Ae^{j\theta} = A\cos(\theta)+jA\sin(\theta) = I + jQ$$
One reason complex numbers (I/Q representation) for RF signals are useful is that for real signals, the positive and negative frequency spectrums are complex conjugate symmetric as complex numbers. This means we can represent any real RF signal using just the positive (or negative) spectrum as an analytic baseband signal, describing the signal independently from the carrier frequency.


Why many DACs uses series of half band filters for interpolation instead of single one? And why, according to some sources, these filters have different number of taps (the last is the shortest one)? Does it really reduce a cost of computation?

The reason cascaded half-band filters, for interpolation applications, should have different numbers of taps can be found at:
https://www.dsprelated.com/showarticle/903.php


Is this a kind of aliasing error and how to make correction? I use a Discrete Fourier Transform to have a representation in the frequency domain:

Notice that, both the real and the imaginary line are softly reaching zero on the left, from 0.5 Hz to 0 Hz. This is a right function.

Subsequently, I wish to multiply the above complex frequency function with another frequency function: DFT*TF, where TF is in the form

The result appears an undesired shift of the imaginary line (and of the module, as a consequence) instead of the previously observed reduction to zero of the above mentioned lines. This shift is a problem because of causing wrong results to the following signal processes.

How could we give a theoretical explanation what kind of problem it is and what have been so far explained about such cases?

We may say that it happens due to the boundary conditions of zero out of the sampling period (the interval I apply the DFT is from 0 to 2047 sampling points, I plot only the part from 0 to 3.5 Hz below) and drives to this shift (see on the left) which is an aliasing error.

Have you ever heard something about such a case?

A practical explanation I could give to this is that, we cannot avoid this phenomenon when multiplying two complex numbers (functions herein) where the imaginary part of (at least) one of them is reduced to zero and is low but not lower than the real part reduction to zero. For example, see below I attach indicative arithmetical values for which the first seven imaginary parts, being a little larger than the real ones, causes a great problem in the following methodology.

The correction I thought is to process the imaginary parts and reduce them to values at least equal to the real ones. Then all the following results works right !
Well, I 've tried this process and works but I cannot be sure myself before asking the experience of others.

What do you think? Could you please indicate me something from books or a useful link?

Complex multiplication dictates that even a complex series with a disappearing imaginary when multiplied with a complex number:
let

$$\begin{align}
z_n &= a_n + i\cdot b_n,\quad a_n, b_n \in \mathbb{R}\\
&\text{then}\\
 z_1 \cdot z_2 &= a_1\cdot a_2+a_1\cdot ib_2+ ib_1 \cdot ia_2 + ib_1\cdot ib_2\\
&=(a_1a_2-b_1b_2) + (a_1b_2+a_2b_1)i\\
&\text{let }b_1\equiv0\\
&=(a_1a_2-b_1b_2) + ({a_1b_2})i
\end{align}$$

and hence, the imaginary part of $z_1$ multiplied with a complex number is not inherently zero, just because the imaginary part of $z_1$ is zero.


What is the $\mathcal{Z}$-transform of a constant? The Fourier transform of a constant exists. Can anyone please tell me what the $\mathcal{Z}$-transform of a constant is? Thanks in advance.

if you define the constant like this:
$$ x[n] = C ~~~,~~~\text{ for all } n $$

Then its $\mathcal{Z}$ transform
$$ X(z) = \sum_{n=-\infty}^{\infty} x[n]z^{-n} = \sum_{n=-\infty}^{\infty} C z^{-n} $$

does not converge for any value of $z$. Hence its $\mathcal{Z}$ transform does not exist.

Note that $\mathcal{Z}$ transform must be an analytic function with continuous derivatives of all orders unlike the DTFT which is not necesarily an analytic function and which would admit the use of impulses (or non-continuous functions) in its expressions.  Then this would allow the DTFT of the constant $x[n] =C$ to be $X(e^{j\omega}) =C 2\pi \delta(\omega)$. But as stated, impulses (which are not analytic functions) are not allowed to exist in $\mathcal{Z}$ transforms.

Had you defined the constant as the practical signal like
$$ y[n] = C u[n]$$
where $u[n] = \begin{cases} 1 &,& \text{ for } n \geq 0 \\ 0 &,& \text{ for } n < 0\\ \end{cases} $ is the unit step function, then its $\mathcal{Z}$-Transform would exist:
$$ Y(z) = \sum_{n=-\infty}^{\infty} y[n]z^{-n} = \sum_{n=0}^{\infty} C z^{-n} = \frac{C}{ 1- z^{-1} } ~~~, ~~~ ROC~~ |z|>1$$

The distinction is clear I suppose.


How to detect curve in chart grid? I would like to detect a curve in a grid. Example (hand-drawn):

https://imgur.com/a/WsE85

I would like to detect the intersects of the curve with the vertical gridlines marked with the "x" symbols. I already managed to get a decent solution when the curve and the grid have different colors (purely based on color separation and an explicit declaration of the horizontal coordinates to measure at).

What I am struggling with is doing the same when grid and curve have the same color (and mostly same color intensity/thickness of lines).

What methods could I look into to solve this? Is there a proven approach for this type of problem?

NB: So far I have worked with OpenCV and Python.

TIA
Matt

Conceptually, I would go for the following steps:

Perform thresholding of the image to create black/white image.
Detect the grid in the image. How you do this depends on the image properties. For example, will the grid lines always be horizontally and vertically or can the grid also be rotated? If it's always horizontal, you can simply add up the pixel values in all rows or in all columns, and you will see the grid lines in the sums. Also, you can detect the grid lines using Hough Transform.
Remove the grid lines from image, e.g. by setting all pixels that belong to the grid to 0.
Use dilation to close the gaps in the curve, where the grid lines were removed.
check, where the curve crosses the grid lines by checking, where the grid image and the dilated image are both 1.


What is the point of decimation? The following makes sense to me:

Decimation reduces the number of samples so it reduces the amount of data that is metered across networks and reduces the processing time for chips. That sounds great.

However, the following does not make sense to me:

I see decimation used in gnuradio-companion all the time.  Often the signal is decimated to reduce the sample rate, then later the sample rate is increased somehow.  What is the point of this?

Sorry for the vagueness of this question.

In communications and in typical GnuRadio applications, decimation is most commonly used to reduce the sampling rate of an oversampled signal, in order to reduce the computational complexity of the system.

Consider this example: You design a stereo FM receiver in GnuRadio. The FM signal is approximately 200 KHz wide, so you sample the downconverted signal at (let's say) 512 kHz. After FM demodulation, the L+R audio signal is present between 0 and 15 kHz; to recover it, you low-pass filter the signal. This leaves you with a signal whose spectrum extends to 256 kHz, but has actual content only up to 15 kHz. It makes sense to decimate it so that its sampling rate becomes, let's say, 32 kHz. In any further processing of this signal, you just reduced the number of operations required by a factor of 512/32=16.

Sometimes you have to interpolate (increase the sampling rate). Let's say that now you're creating a stereo FM signal. You need to shift the spectrum of the L-M signal from the range 0-15 kHz to the range 23-38 kHz. If the signal is originally sampled at 32 kHz, the operation is impossible. You need to first interpolate it to a new sampling rate of at least 38*2=76 kHz.

(See https://en.wikipedia.org/wiki/FM_broadcasting#Stereo_FM if you're unfamiliar with the format of FM broadcast signals).


Count objects in the image How an algorithm can be developed to count numbers and bolts in the image below.

I will be using MATLAB as part of my solution.  The basic algorithm will be this:

Read in the image and convert to black and white
Invert the intensities and fill in the holes.  There are speckles in the image where when you convert the image to black and white, the objects are not solid.  We want to ensure these are solid
Count how many unique objects there are using bwlabel.

You'll need the Image Processing Toolbox for this algorithm.  If you don't have this and are using MATLAB, let me know and I'll edit my post.

Without further ado:

% Read in image, convert to black and white - Link comes from your image posted here
im = imread('http://i.stack.imgur.com/lBGU1.png');
imBW = im2bw(im, 0.3); %// Specify manual threshold of 0.3

% Invert intensities and fill in holes
imBWFilled = imfill(~imBW, 'holes');

% Count how many unique objects there are
[L,num] = bwlabel(imBWFilled);

% Show final image and display number of objects counted in the title
imshow(imBWFilled);
title(['Total number of objects: ' num2str(num)]);

L contains a map where each pixel is an ID of which object that belongs to.  0 means that the pixel belongs to the background while any pixel that has values greater than or equal to 1 means that the pixel belongs to a particular object associated with that ID number.  num gives you the total number of objects seen in the image.  The output of this code thus gives:

If you don't have MATLAB and want to compute this using another language, let me know.  You can find the number of objects by considering an image as a connected graph.  If you don't have bwlabel, you can use any graph searching algorithm (breadth-first search, depth-first search, etc.) to help you compute this.  Start with any pixel that belongs to an object, and perform BFS / DFS to visit those pixels that belong to the object.  You then set all of these pixels to belong to an ID number.  When the queue / stack is empty, you then select another pixel that belongs to an object and repeat the algorithm.  You stop when you have visited all pixels that belong to objects.  The total number of IDs you have issued will essentially be how many unique objects you have in your image.

Edit - July 14th, 2014

Seeing your comments, you wish to be able to distinguish between what is a nut and what is a bolt.  As such, we can simply add on top of the current code.  The algorithm that I have developed is loosely based on circularity:

For each object, find the centre of mass.  This is simply taking all of the X and Y co-ordinates and averaging them.
For each object, find their boundaries / perimeters
For each centre of mass we have, calculate the distance between this point and all of the object boundary points.
Find the difference between the maximum and minimum of these distances - denote this the range
We will have N ranges for N objects.  Simply take a look at these ranges and threshold the array.  If you take a look at the bolts, they have a longer vertical direction than horizontal direction.  For the bolts, the vertical and horizontal direction is pretty much the same.  As such, the range should give us an idea of what is a nut and bolt because the difference between the maximum distance to the centre and the minimum distances should be very large in comparison to the nut.  As such, for each object we have, check to see if the range is below a certain value (nut) or above a certain value (bolt).

Without further ado, here is the code:

%%% For each object in the image, find the centre of mass
centres = zeros(num,2);
% Cycle through each unique object label and extract (X,Y) co-ordinates
% that belong to each object.  Compute centre of mass for each.
for n = 1 : num
    bmap = L == n;
    [rows,cols] = find(bmap == 1);
    centres(n,:) = [mean(cols) mean(rows)];
end

% Find boundaries of all objects
bwBound = bwperim(imBWFilled, 8);

% For each object, find the distances between the centre of mass with all
% of the pixels along the boundary for each object.  Find the range (max -
% min).
ranges = zeros(num,1);
for n = 1 : num
    bmap = L == n; % Obtain all pixels for an object
    boundPix = bwBound & bmap; % Logical AND with boundaries map to extract
                               % only those pixels around the perimeter
    [rows,cols] = find(boundPix == 1); % Find these locations
    % Compute the distances between the centre of mass with these points
    dists = sqrt((cols - centres(n,1)).^2 + (rows - centres(n,2)).^2);

    % Find the difference between the maximum and minimum distances
    ranges(n) = max(dists(:)) - min(dists(:));
end

This isn't complete yet.  I stopped here so you can see what the ranges for all of the objects look like:

ranges =

37.9615
63.9613
54.9266
 5.0716
 4.1578
 6.3114
 7.2356
41.6381
10.3123
34.0938
 5.0021
67.3290

As you can see, there are some ranges that are quite small ($< 10$) while there are some ranges that are large.  To be safe, let's choose a threshold of 15.  As such, for those ranges that are less than 15, this is classified as a nut while those that are larger are classified as a bolt.  When we're done, we simply count how many of these fall within each range, and those are the number of nuts and bolts we have.  Let's continue with the algorithm:

% Find those object IDs that have less than a range of 15.  These are the
% bolts
% The rest are nuts
indBolts = find(ranges < 15);
indNuts =  find(ranges >= 15);

% Total number of nuts and bolts
numBolts = numel(indBolts);
numNuts = numel(indNuts);

This will give you 6 and 6 respectively as you expected.  Now to finish up the algorithm, I'm going to give you an added bonus.  For each unique object we have, I'm going to colour their interiors to be a certain shade of gray.  The background will be black, while a nut is gray and a bolt is white.  The code to do this is:

finalMap = uint8(zeros(size(imBW)));

for n = 1 : numBolts
    finalMap(L == indBolts(n)) = 128;
end
for n = 1 : numNuts
    finalMap(L == indNuts(n)) = 255;
end

figure;
imshow(finalMap);
title(['Number of Nuts: ' num2str(numNuts) ', Number of Bolts: ' num2str(numBolts)]);

The final image we get is:


Octave error - STFT parameters when detecting pitch in guitar audio signals I am working with this STFT method to detect pitch in monophonic guitar audio signals.

The method is working correctly for medium pitches (above 150Hz). However, in low notes I have some issues.

For example, I have an audio file that had the A2, A#2 and B2 notes played. So, the pitches detected should be approximately 110, 116, 123. When I ran the STFT with 40000 as a sampling rate, there was an octave error in all of them. The notes detected were A3, A#3 and B3.

I though I should see why this happens and I went in my STFT, and printed, for each note played, the 5 strongest bins in terms of magnitude. Those were the results:

First note: [539.49506, 975.88092, 419.11325, 326.00168, 217.16208]

Second note: [0.0, 0.0, 573.76556, 229.9518, 117.29102]

Third note: [0.0, 240.0067, 0.0, 118.64342, 0.0]

As you can see, for the first note, the correct pitch is not even amongst the 5 max. magnitudes. For the second and third notes, they are amongst the top 3 but not enough to become the selected pitches.

I tried to do the same with a sampling rate of 30000, and I still had the same problem, but this time the correct pitch was in the top 5 for the first note.

I know how to fix this problem, so my question is more theoretical. How can the sampling rate affect this? Why did this happen with 40000 and not 30000? In the end, how should I decide on my SR final value?

As a follow-up question, I see that I can also choose on other parameters like the n_fft (number of FFT bins or FFT size) and hop_length. How should I choose on those? What are upsides and downsides for high/low values (if the answer is too long to answer, then a reference would be great).

The octave error is because you are using a spectral frequency estimator (STFT), not a psychoacoustic pitch detection or estimation method.  The spectral frequency peaks in the sound are often not the pitch frequency, but instead often one of the many harmonics.  As you found, the fundamental pitch frequency can even be completely missing.  This is especially true for guitar notes, voices, and sounds from other large stringed instruments.

The solution is not to play with sample rate and STFT FFT window parameters, but to look up more reliable pitch detection and estimation algorithms.


At which point in the signal processing chain to use Hann Window? I have a question about signal processing. Experienced with embedded design and code, but DSP is in new territory.

Thus far I am using the CMSIS DSP library functions, which are working well. I am sampling a low frequency signal. Using a timer interrupt, the sampling rate is 40x higher than the highest signal I expect to see. The "buffer" size I process is 2048. I first apply a low pass FIR (arm_fir_f32), and then to the arm_cfft_f32 function so I can then calculate the frequency information. It works well enough to a point where the frequency is +/- 1Hz. However, I am hoping to get this to +/-200mHz.

I understand that I need to reduce the power (noise) at the endpoints, and the Hann window should be the best (better?) choice. Thus my questions:

At which point do I apply the Hann algorithm? (as follows)

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

Should $N$ be the size of my sample block?

EDIT:

Code for windowing function:

(INPUT_SAMPLES=2048)
for(uint16_t i = 0; i < INPUT_SAMPLES; i++)
{
    ADC_Values[i] = 0.5 - (0.5 * cos ( (2.0 * PI * ADC_Values[i]) / (INPUT_SAMPLES - 1)));
}

Thank you for any help on this,
Roger

EDIT 2:

jojek,

creating the window and changing the calculation as suggested has fixed it. Although the window has not improved resolution as much as I'd like, it's better, and close enough to be quite usable. Fantastic...thanks for all the replies everyone.

You should apply window in time domain. The purpose is to make your signal go to zero at boundaries of your frame/sample block. It is necessary because DFT assumes signal to be periodic, and zeroing it at the ends makes it so. Obviously if length of your frame is constant, then precompute window values.

Enough words:

It just occured to me that you want to get the frequency resolution of 0.2 Hz - is that correct? You must keep in mind that frequency resolution is given by: $\Delta f = \frac{f_s}{N} $. So unless your sampling frequency $f_s$ is not lower than approx. 410 Hz (unlikely) you won't achieve the resolution of 200 mHz. You need to extend your observation window. Also padding your signal with zeros will not help too much - this will be only an interpolation of your spectrum points.

EDIT:

If you perform any FIR filtering, then apply the window afterwards. Windowing should be done just before the FFT calculation to make your signal 'disappear' at the ends (especially when filtering produces artifacts for first output samples). Regarding your frequency resolution, with 2kHz of $f_s$ and $N$ equal to 2048 you get approximately 0.98 Hz of a resolution. If  you want more than that, then you might think about extending your sample block or decreasing the sampling frequency. Just keep in mind equation I mentioned above.

EDIT 2:

I presume that the function you pasted (and one I moved to your original question) is only for calculation of window. Do I understand you correctly? Or is it function that is performing windowing of your signal (variable ADC_Values[i] is present). Check the following things:

I corrected your equation for window function - notice that it is
dependent only on n - there are no signal samples.
If you want to perform windowing within your function then code
should look more like:

ADC_Values[i] = (0.5 - (0.5 * cos ( (2.0 * PI * i) / (INPUT_SAMPLES - 1)))*ADC_Values[i];

Window calculation itself is given by:

`win[i] = (0.5 - (0.5 * cos ( (2.0 * PI * i) / (INPUT_SAMPLES - 1)));`

Therefore I suggest you to precompute the win once and for all, passing it later to your function and performing ADC_Values[i] *= win[i];. This will save some time.

I am assuming that you making typecasting correctly. If you have any doubts, then extract your signal block before and after filtering, together with window samples and windowed signal. This should help you to debug the code.

EDIT 3:

I am afraid that you still do not understand the principle. Changing of window can only remove leakage effects - kind of noise covering your harmonics. It will not increase your frequency resolution. In the matter of fact, it will make it worse - as rectangular window has narrowest main-lobe. Let me repeat once again - if you want to increase your frequency resolution then you must use the appropriate buffer size and sampling frequency, according to equation $\Delta f = \frac{f_s}{N} $.


Deriving 2-D discrete Fourier transforms I have a problem in DFT. It was one of my past-year exam papers questions.

Question:

  Let $F(u,v)$ be the 2-D Fourier transform of a 2-D continuous function
  $f(x,y)$. Derive in terms of $F(:,:)$ the 2-D Fourier transform of
  each of the following functions

  1) $f(x,-2y)$

  2) $f(x+2y,y)$

I know how to do 1-D Fourier transforms but not 2-D. I'm not sure how to start on it and need some guidance.

For the second part, this was my approach. Please let me know if it's right or correct me if it's wrong.

Let $τ= x + 2y$ hence $x = τ-2y$ and $dx = dτ$
$$
\begin{align}
\mathfrak{F}\{f(x+2y,y)\}&=∬ f(τ,y)e^{−j2π(u(τ-2y) +vy )} dx\ dv\\

\mathfrak{F}\{f(x+2y,y)\}&= ∬ f(τ,y)e^{−j2π(uτ+(-2u+v)y  )} dx\ dτ \\
\mathfrak{F}\{f(x+2y,y)\}&= F(u,-2u+v)
\end{align}
$$

Here is the first one:

By definition,

$\mathfrak{F}\{ f(x,-2y)\} = \iint_{- \infty}^{\infty} f(x,-2y)e^{-j2\pi (ux+vy)}dxdy$

let $\tau = -2y$ and conversely $y=\frac{-\tau}{2}$

$\mathfrak{F}\{ f(x,-2y)\} = \iint_{- \infty}^{\infty} f(x,\tau)e^{-j2\pi (ux-\frac{v\tau}{2})}dxd(-\frac{\tau}{2})$

$\mathfrak{F}\{ f(x,-2y)\} = -\frac{1}{2}\iint_{- \infty}^{\infty} f(x,\tau)e^{-j2\pi (ux-\frac{v\tau}{2})}dxd\tau$

$\mathfrak{F}\{ f(x,-2y)\} = -\frac{1}{2} F(u,-\frac{v}{2})$

The other one will be harder, but I will leave it to you.


Theoretical Maximum of DFT For any discrete input signal between +1 and -1, what is the theoretical maximum DFT?

If the input is a cosine of $N$ samples with amplitude $A$, the peak spectral magnitude is $A*N/2$. But what about any arbitrary waveform?

A square wave with the same amplitude as the cosine is composed of various sines, the fundamental frequency of which has a greater magnitude than the square wave. Therefore the Fourier transform of a square wave with amplitude also between $1$ and $-1$ is greater than the cosine with amplitude between $+1$ and $-1$.

But generalizing this to any waveform, ie a triangle wave, sawtooth wave, and also non geometrically "nice" waves, all with amplitudes between $+1$ and $-1$, what is the theoretical maximum?

I want to know because I need to represent the output of a DFT with 16 signed bits. However if the magnitude of the DFT is greater than $2^{15}$ I would need to scale. I therefore need to determine if this will ever happen (for any input waveform) and how much the scaling factor should be.

Thank you very much,

A simple bound is just  the number of samples. Indeed, for
$$ X_k = \sum_{n=0}^{N-1} x[n]e^{-2i\pi kn/N}$$
then obviously:

$$ |X_k| \le \sum_{n=0}^{N-1} \left|x[n]\right|\left|e^{-2i\pi kn/N}\right|$$
hence
$$ |X_k| \le \sum_{n=0}^{N-1} \left|x[n]\right|\le N$$
if $-1\le x[n]\le 1$. The $N$ bound is
s attained for a given $k$ when you choose a complex $x[n]=\overline{e^{-2i\pi kn/N}} = e^{2i\pi kn/N}$. If your input allows a constant signal  (equal to one), $x[n]=1$, then this is attained for the zeroth frequency $X_0 =  \sum_{n=0}^{N-1} 1 = N$. As a side note, if the signal length is a power of $2$, this is quite convenient.

For more restricted inputs, like real, zero-average signals only, there could be a tighter bound. A recent paper presented at SampTA 2017 provides DFT bounds for some random sequences (Bernoulli): Bounds on Discrete Fourier Transform of Random Mask (published version).


how a signal for which calculate RMS can be filtered and have a fast settling time? i designed a sinusoidal pwm amplifier (3KHz to 20KHz) that regulates its output voltage reading the RMS value of output current.

The RMS is obtained squaring the signal sampled by AD converter @200KHz and then applyng a FIR filter with 120 TAPS.

The advantage of this method is that the settling time of the FIR (so the response to amplitude changing) is about 600us so i can correct the pwm duty every 1ms because the value is alread stable.

All works fine but the problem arises when an external signal is added to the signal the amplifier is generating.

Infact a signal of about 600Hz different from original one (i.e. amplifier generates 4950Hz and a 5550Hz signal is added as external noise) causes an oscllation in RMS measure makes the amplidifer works in wrong way.

Applying a bandpass FIR makes the response of FIR slower. I tried downsampling by a factor of 4 and make a FIR of 1KHz bandwith but the fastest settling time i obtained was about 5ms that is too much because i've to regulate every 1ms.

Have you any idea or advice? Is there any algorithm to cancel unwanted component maintaing fast RMS response?

Thank you . Let me know.

I believe you are running into stability issues in an output power control loop design.
See below a diagram of similar power control loops that I have implemented, where for stability reasons any filtering in the loop is minimized and only done with the loop filter itself which is designed for stability. The noise you are trying to filter gets filtered by the loop bandwidth itself.
The diagram below is very simple to implement (first order type 1 loop), where the accumulation of the derived error sets the loop bandwidth by adjusting K.  Make K very small for a long averaging time and better filtering but slower response.

Ultimately if the trade of filtering versus response time with a simple implementation like shown above is insufficient, you can change the loop filter to be a PID loop. This would more complicated and require more understanding of control loop design and a more detailed loop model but give you better performance in noise filtering versus response time. Essentially the accumulator (I) shown above would also be summed with a proportional path (P) and a derivative path (D) all with different gain weights and followed by a second accumulator for a second order type 2 PID loop.
Once a loop model is created, the modeling of the response time and effective filtering is very straightforward with tools such as Matlab, Octave or Python. As an example of the design process, below shows the loop model for the Power Control Loop given by the block diagram above along with the modeled response times versus K:

The above case results in an open loop gain given by:
$$G_{OL}(z) = \frac{0.83155K}{z(z-1)}$$
Below demonstrates given the Matlab and Octave code snippets how simple it is once the model is established to create the step response for the loop. Note the overshoot indicative of a higher order loop- this is due to the parasitic $z^{-1}$ in the digital loop implementation, and with a high enough K the loop will be unstable (unlike a true first order type 1 loop). Similarly the noise bandwidth, etc is easily established and optimized.


DAC Sinc Roll-off vs Freq The output of a Digital to Analog converter has a magnitude roll-off in frequency as a Sinc function with the first null at the sampling rate. Why is this the case?

Not ALL DACs (D/A converter) do this, but a traditional DAC operates on the digital samples as a zero-order hold; each sample is held in the analog at a constant value until the next sample. This process is identical mathematically to convolving the discrete samples (as weighted impulses) with a rectangular function of width T where T is the sample time.

The Fourier Transform of a rectangular function of width T is a Sinc function with first null at 1/T, where 1/T is the sampling rate.

Convolving in the time domain is equivalent to multiplying in the frequency domain, so the spectrum of your discrete signal is multiplied by the Sinc function in the process of analog reconstruction in the DAC due to the zero order hold.

There are other DAC implementations that reduce this effect such as interpolating DACs and "Return to Zero" (R2Z) DACs. An R2Z DAC forces the output to return to zero prior to the next sample, reducing T and therefore pushing out in frequency where the first null occurs (therefore reducing passband rolloff).


How do Phons relate to loudness? Let's say I have the very same signal stored in 2 .wav files. And lets say I process each file using a function, and it is established that one file is 20 phons larger than the other. Lets ignore how that function calculates this.

Does this mean that one file is approximately 4 times louder than the other?

What formula can I use to relate the difference in phons for two signal to the difference in their loudness?

-10 phons = 0.5
0 phons = 1
10 phons = 2

Your equations are correct (as a first order approximation).

The phon is an attempt to establish a perceptual loudness measurement. Perceived loudness is rather complicated since it needs to take into account what happens in the human auditory system and in the human brain.

Here is is how it works: First you measure the sound energy at the listener's position. The relevant physical quantity is the variation in air pressure, called "sound pressure" and measured in Pascal or Pa like any other pressure. This is then referenced to 20 micro Pascal and turned into dB like this
$$ x = 20\cdot log_{10}\frac{p}{20\mu Pa} $$
where p is the measured sound pressure and x is the sound pressure level in dBSPL. This is a direct measurement of the physical sound energy.

The human auditory system has different sensitivity at different frequencies. It's most sensitive around 2kHz-4kHz, a lot less sensitive at 100 Hz and at 30 kHz you don't hear anything at all, regardless of how much sound energy is there. The "phon" takes this into account. It's based on measurements that compare the perceived loudness of sine waves at different frequencies. The results are often referred to as "equal loudness contours" (see for example http://en.wikipedia.org/wiki/Equal-loudness_contour).

A loudness of X phons means "as loud as a 1 kHz sine wave at X dBSPL". For a 1 kHz sine waves phons and dBSPL are the same thing. For sine waves of other frequencies phons and dBSPL are related through the equal loudness contours. For more complicated signals calculating phons is a bit more tricky.

Phon is not a direct loudness measurement. Perceived loudness, as most human modalities, follows roughly the Weber-Fechner law (e.g. http://en.wikipedia.org/wiki/Weber%E2%80%93Fechner_law). As a first order approximation it takes ten times the energy to double the perceived loudness (all other things being equal). An increase in 10 phone doubles the loudness and an increase in 20 phone does, indeed, increase the perceived loudness by a factor of 4.

To calculate phon (relative or absolute) from a wav file, the file MUST be calibrated in units of Pascal (sound pressure) since the equal loudness contours are non linear and change with absolute level.


Hilbert transform with scipy.signal I'm coding with python and I just started to use the Hilbert transform from scipy.signal. Consequently, to check the result I got I tried to verify a property of the Hilbert transform:

$H(H(\vec{x})) = -\vec{x}$

So I implemented the following code:

import numpy as np
from scipy.signal import hilbert

# my input vector, an 1D numpy array of real numbers
x = ...

# I consider the imaginary part here as scipy.hilbert
# returns the analytical signal computed with the
# Hilbert transform (imaginary part).
x_ = hilbert(hilbert(x).imag).imag

However I can't retrieve x_ == -x.

Do you know what this could be due to?

Here are plots of x and x_ (first figure) and the difference (second figure):

What you did should work, up to a constant, as mentioned by Dan Boschen in a comment. The Hilbert transform will remove any DC component, so the output of the second Hilbert transform should equal the (negative of the) input, up to a constant, and, of course, up to numerical accuracy.

The following example (in Octave) shows this:

t = 0:1024;
x = sin(pi/7*t)+3*cos(pi/5*t);    % some signal
xh = imag( hilbert(x) );
xhh = imag( hilbert(xh) );
err = x - mean(x) + xhh;
max( abs(err) )                   % 6.2172e-15


Projecting a vector to another to detect the sparse values of such vector Assuming we have sparse vector of length $N$ such as $X = [0,1,0,-1,0,1,1,0]$ which has some non-zeros values. The vector $x = iFFT(X)$ is convoluted with another vector $h$ resluting $y = h*x$. Suppose that $h$ and $y$ are known, is it possible to recover the sparse vector $X$ using such methods of sparse vector estimation, such as compressive sensing or any other method?

$NP$:
I have suggested to take $FFT$ for $y$ resulting $Y = FFT(y) = HoX$, where $o$ is the element-wise multiplication between $H$ and $X$. Then based on $Y$ I think it's easier to recover sparse vector $X$. But the issue when using element-wise division, I can't get the position of each non-zeros element in $X$. On the other hand, I don't know how to implement compressive sensing to recover $X$.

Thank you in advance

When you have the described convolution measurement, you can indeed recover your $X$ using sparse recovery methods. Some examples are, as you suggest, OMP, MP, or for example subspace pursuit (SP), compressive sampling matching pursuit (CoSaMP), iterative hard thresholding (IHT), a plethora of other variants of greedy recovery methods, as well as basis pursuit (BP), basis pursuit denoising (BPDN), and approximate message passing (AMP).
The key to this is the measurement matrix. Have a look at the random demodulator architecture to see how they construct their measurement matrix ($\mathbf M$) [1]. This is very similar to what you want.
Your IFFT is modelled by their $\mathbf F$; in their case the transform is just DFT instead of IDFT, but that does not change the principle. Your measurement matrix (lets us call it $\mathbf H$) is a little bit different than [1]'s $\mathbf M$.
Assume that we can express your convolution as discrete convolution. I would model you measurement matrix as:
$$\mathbf H = \begin{bmatrix} \mathbf h & 0 & & \ldots && 0\\
0 & \mathbf h & 0 & \ldots && 0\\
\vdots\\
0 & & \ldots && 0 & \mathbf h
\end{bmatrix}$$
This is a convolution matrix. The non-zeros parts $\mathbf h$ of each row of $\mathbf H$ is the vector you convolve $x$ by (where I would call it vector $\mathbf x$). The non-zero parts of $\mathbf H$ generally overlap across the columns.
See for example MATLAB's convmtx for a practical example of construction of $\mathbf H$. You may want to sub-sample the rows of $\mathbf H$ (i.e. only keep every $m$th row), corresponding to a lower sample-rate of your sequence $y$ than your original sequence $x$.


Frequency Domain Distribution I have a complex signal in the time domain normally distributed. What will be its distribution in the frequency domain?
I assumed since the frequency domain is a linear transformation the distribution will not change.

I assume you mean that at each time $t$, the signal is a normally-distributed random variable.  This tells you the probability that the signal will take a value in any given range, but it does not tell you whether/how the signal values at different times are related to one another.  Typically, for a random signal one defines not only the probability distribution for the signal at each time, but also the autocorrelation function $E \{ x(t) x(u) \}$ or some other measure of how signal values at different times are related or independent.  This is required in order to determine how the signal is distributed in the frequency domain.


Transfer function sinusoidal response Why does the phase shift between the input and the output of a transfer function vary with the frequency of the input sinusoid?

For any LTI (Linear Time Invariant) system (considering the steady state response i.e. ignoring the transient response which generally becomes negligible in a very short period), if the input is a sinusoidal signal, the output is always a sinusoidal signal of the same frequency.  But the amplitude of the output sinusoidal can be different from the input amplitude.  And the ratio of the output amplitude to input amplitude is a constant for that particular frequency - whatever the input amplitude and phase are. Similarly, the difference between the input and output sinusoidal phases is also constant for the particular frequency - irrespective of the amplitude. The transfer function can be used to give this input-output amplitudes ratio and input-output phase difference as a function of frequency.  Hope this helps.  The frequency dependency is derivable from the input/output differential equation in analog case or difference equation in discrete case.


How FIR filters provide a linear phase? I know that we can design a linear phase FIR filter by choosing filter coefficients symmetric or antisymmetric, what is the intuitive idea behind that?

The reason is Euler's formula, from which you get

$$\cos(\omega)=\frac12\big(e^{j\omega}+e^{-j\omega}\big)\tag{1}$$

and

$$j\sin(x)=\frac12\big(e^{j\omega}-e^{-j\omega}\big)\tag{2}$$

If you have symmetric or anti-symmetric coefficients, the corresponding frequency response can always be decomposed in purely real-valued cosine terms $(1)$ or purely imaginary sine components $(2)$, plus a linear phase term if the filter is not centered around $n=0$.

Take as a simple example a filter of length $N=5$ with filter coefficients

$$\mathbf{h}=\big[1, 2, 3, 2, 1]\tag{3}$$

starting at index $n=0$. The corresponding frequency response is

$$\begin{align}H(e^{j\omega})&=1+2e^{-j\omega}+3e^{-2j\omega}+2e^{-3j\omega}+e^{-4j\omega}\\&=e^{-2j\omega}\big[e^{2j\omega}+2e^{j\omega}+3+2e^{-j\omega}+e^{-2j\omega}\big]\\&=2e^{-2j\omega}\big[\cos(2\omega)+2\cos(\omega)+3\big]\tag{4}\end{align}$$

Since the last term in brackets is purely real-valued, $H(e^{j\omega})$ has a linear phase response $\phi(\omega)=-2\omega$.

For even-length filters you get an additional delay of half a sample, and for anti-symmetric filters there's an additional phase shift of $\pi/2$ due to $j$ on the left side of $(2)$.

Also take a look at this answer discussing the four types of linear phase FIR filters.


Designing an FIR filter in Matlab with differing number of precursor taps and postcursor taps I'm trying to fit an FIR filter to the measured frequency response (magnitude and phase) of a device. The frequency response is essentially a low-pass filter. I'd like to vary the number of taps to see how many are necessary to appropriately match the response, but also give differing number of precursor and postcursor taps in the design.

I've investigated Matlab's fdesign.arbmag and fdesign.arbmagnphase functions, but these seem to only give designs with the same number of pre- and post-cursor taps.

Here's some sample code with simplified frequency (F) and amplitude (A) values for reference.

F = 0:0.1:1;
A = [1 0.85 0.9 0.7 0.5 0.3 0.35 0.1 0.2 0.1 0];
filterOrder = 7;

d = fdesign.arbmag('N,F,A',filterOrder,F,A);
Hd = design(d);
fvtool(Hd)
tapWeights = [Hd.numerator];
disp(tapWeights)

The above code runs happily, and produces a low-pass filter that matches the frequency response reasonably well, with tap weights as follows:

[0.0030  0.0361  0.2161  0.4500  0.2161  0.0361  0.0030].

However, the central tap (the cursor, 0.45) is surrounded by three precursor and three postcursor taps. Is there a way to design, say, with one precursor and six postcursors.

There is no reasonable way to explicitly choose the number of taps before and after the peak. The reason is simple: the arbitrary magnitude response design results in a linear phase, and, consequently, the impulse response is symmetrical. If you specify a desired complex frequency response in terms of magnitude and phase, the location of the peak is an outcome of the design process. The smaller the average delay (implicitly specified by the desired phase response), the further to the left the peak will occur.

Take as an example the design of an FIR low pass filter with $21$ taps. The figure below shows $3$ designs. The difference is the desired phase response. The top figure shows the impulse response (filter taps) for a linear desired phase response, the middle one for a non-linear phase response with a lower pass band delay than the linear phase response, and the bottom figure shows the impulse response of the corresponding minimum phase filter. The location of the peak is entirely determined by the specified phase response.


Overshoot computations not matching Take the following transfer function of a 3rd order system:

$$H(s)=\dfrac{2.302~s+0.3548}{s^3+0.739~s^2+3.223~s+0.3548}$$

with poles:

         Pole              Damping       Frequency      Time Constant
                                       (rad/seconds)      (seconds)

 -1.13e-01                 1.00e+00       1.13e-01         8.89e+00
 -3.13e-01 + 1.75e+00i     1.76e-01       1.78e+00         3.19e+00
 -3.13e-01 - 1.75e+00i     1.76e-01       1.78e+00         3.19e+00

and with the following unitary step response:

If I compute the percentage overshoot (PO) based on the damping ratio $\zeta=0.176$, I get:

$$PO=100{\times}e^{\dfrac{-\zeta\pi}{\sqrt{1-\zeta^2}}}=100{\times}e^{\dfrac{-0.176\pi}{\sqrt{1-0.176^2}}}=\boxed{57.02\%}$$

However, if I compute the PO using the graphical method (comparing the peak value with the final value) I get a completely different result:

$$PO=\dfrac{v_{peak}-v_{final}}{v_{final}}{\times}100=\dfrac{1.2-1}{1}{\times}100=\boxed{20\%}$$

I don't understand why such discrepancy. Why are my PO computations not matching?

Computing the partial fraction expansion,

$$H (s) = \frac{2.302 s + 0.3548}{s^{3} + 0.739 s^{2} + 3.223 s + 0.3548} \approx \frac{2.2861 - 0.0309306 s}{s^2 + 0.626454 s + 3.1525} + \frac{0.0309306}{s + 0.112546}$$

For the time being, let us neglect the step response of the 1st order subsystem. Note that the 2nd order subsystem has the following form

$$\pm \gamma (s - z) \left(\frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}\right)$$

where $\gamma \in \mathbb R$ is the gain (or attenuation) and $z \in \mathbb R$ is a (finite) zero. However, the formula you used to calculate the step response overshoot is only applicable for 2nd order systems of the form

$$G (s) := \frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}$$

whose zeros are at infinity. To summarize, the overshoot is smaller than you expected because:

you neglected the 1st order subsystem.
you neglected the attenuation $\gamma \neq 1$.
you did not subtract $\gamma$ times the derivative of the step response of $G (s)$ at the peak time.

In any case, not neglecting the above would only provide one estimate of the overshoot. To compute a satisfactory approximation of its exact value, you can take the inverse Laplace transform, differentiate and find where the derivative vanishes.


Find beginning of signal for BPSK Demodulation I want to store a small amount of data in a sound file recorded by a microphone but I know nearly nothing about signal processing. I managed to demodulate my data encoded using binary phase shift keying with this schema:

(From here: http://courses.cs.washington.edu/courses/cse561/10sp/project_files/FinalReport_FCI_YH.pdf)

To find the beginning of the signal I tried searching for the time where the strength of my carrier frequency in the signal exceeds a certain threshold but this is way too unreliable.

What would be a good way to solve this?

This problem is typically solved by preceding your actual payload data with a synchronization pattern, or sync word. This is a sequence of symbols that are known to the receiver ahead of time and are used to aid synchronization with the payload data frame. It could be detected in a couple ways:

Sync pattern matched filtering: When looking for the synchronization marker, instead of just integrating over one data bit, you could create a matched filter that represents the entire modulated sync word. This has the advantage that the SNR at the output of the sequence matched filter is higher than you would have otherwise, by a factor of the number of bits in the pattern. This can help if you're trying to detect the frame start in a very low-SNR environment.

Disadvantages of this approach are that it is more computationally complex (because you're correlating against an entire sequence of bits), and also more sensitive to frequency offset between the transmitter and receiver.
Hard decision pattern matching: Instead of the "soft-decoding" approach of correlating against the entire modulated sync sequence, you could instead just look for the pattern's sequence of symbols at the output of your existing receiver. As an example, you could decide that all of your frames will begin with the sequence 01001100101010 (I make no warranty that this is a good choice). You would then slide through the bits that come out of your BPSK decoder and declare the start of a frame when you see the sync pattern. It may also help to allow some small number of bit errors in the pattern when looking for the sequence.

Compared to coherent integration across the entire sync sequence, this technique will not perform as well in low-SNR environments. However, it is much simpler and may work just fine for your application.


How do you program a completely transparent crossover? Or, if that's too broad, what is/are the most popular algorithms?

Background: I have no formal DSP training but much informal tinkering. I am trying to program a crossover for an audio effect. Another audio effect I've used (ReaXComp) has completely transparent crossovers, or so it sounds to my ears. I assumed that crossovers were just two matching filters, a lowpass and highpass filter of some kind, at the same frequency with a Q value of $\sqrt{\frac{1}{2}}$ (that seems to be the "default" Q value). However, when I run an audio signal through a lot of such crossovers it sounds... weird.

Later, I discovered through experimentation that you can put a two-pole low pass filter on a sound, and then subtract the original sound from the output of that and it sounds like a high pass filter, which is exatly what I want: a low and high pass filter that when added together create mathematically the exact same sound. Unfortunately, if you want a steeper filter and apply an identical filter again to get a 4-pole filter (again at a q value of $\sqrt{\frac{1}{2}}$) the resulting high band just starts to sound weird, not like a high pass filter at all. It basically sounds like the original sound with some sort of weird comb filtering. So I finally decided to come here and ask how to do it properly.

Edit:

robert bristow-johnson asks in the comments:

okay, RealXComp is a multiband compressor. so you want to split the audio into subbands for multiband compression and, if none of the individual-band compressors hit their knee, you want it all to add up to the original signal (perhaps with a little delay)? is that it? is that what you're trying to do?

Yes, Robert, that's exactly what I'm trying to do, sorry if that wasn't clear, I'm a bit of a DSP noob

I'm no expert but I'll explain how I would approach this. In general, crossover design involves selecting your respective cutoff frequency (or frequencies, in your case) and desired attenuation or filter order, and deciding on a filter type. Trivially, a Butterworth filter is a popular choice as they are easy to design and implement. According to Wikipedia, a better choice may be the Linkwitz-Riley filter:

  A Linkwitz–Riley (L-R) filter is an infinite impulse response filter used in Linkwitz–Riley audio crossovers...It is also known as a Butterworth squared filter. A Linkwitz-Riley "L-R" crossover consists of a parallel combination of a low-pass and a high-pass L-R filter. The filters are usually designed by cascading two Butterworth filters, each of which has −3 dB gain at the cut-off frequency. The resulting Linkwitz–Riley filter has a −6 dB gain at the cutoff frequency. This means that summing the low-pass and high-pass outputs, the gain at the crossover frequency will be 0 dB, so the crossover behaves like an all-pass filter, having a flat amplitude response with a smoothly changing phase response. This is the biggest advantage of L-R crossovers compared to Butterworth crossovers, whose summed output has a +3 dB peak around the crossover frequency.

Wikipedia on L-R Filters

I would design these such that the bands of interest are all separated out nicely, either through some sort of series configuration or by explicitly designing bandpass filters instead of lowpass/highpass filters (except for your lowest and highest bands of interest, of course).

This could either be done by designing the filters as analog and then using a bilinear transform to bring them into the Z-domain, resulting in a form that could be implemented in your language of choice, OR could be done digitally from the get-go in Matlab or Python (using SciPy's iirdesign() function).

SciPy's IIR Design Function

Then, you would need to implement them in whatever language/platform you're using (although I would hazard a guess that its C++ on Windows/Macfor a VST?).

Additionally, the concept of a "transparent" filter is something of a oxymoron. Either you are going to affect your signal or you aren't, and everything will involve some compromise. For instance, better filtering in the frequency domain may result in poorer performance in the time domain. Choosing what you can "ignore" is key to a good design.


What is the difference between neural network and deep neural network. an example thank you What is the difference between neural network and deep neural network in speaker recognition .Can i have  an example of code sur Matlab thank you

The difference resides in the number of perceptron neuron layers a nettwork is build on.

Neural Networks have as its base a perceptron. This is like a human neuron. The perceptron can weigh up different kinds of evidence in order to make decisions. And it should seem plausible that a complex network of perceptrons could make quite subtle decisions.

Neural networks have different layers of perceptrons. First layer of perceptrons make very simple decisions, by weighing the input evidence. What about the perceptrons in a second layer? Each of those perceptrons is making a decision by weighing up the results from the first layer of decision-making. In this way a perceptron in the second layer can make a decision at a more complex and more abstract level than perceptrons in the first layer. And even more complex decisions can be made by the perceptron in the third layer. In this way, a many-layer network of perceptrons can engage in sophisticated decision making.

In fact, a small change in the weights or bias of any single perceptron in the network can sometimes cause the output of that perceptron to completely flip, say from 0 to 1. That flip may then cause the behaviour of the rest of the network to completely change in some very complicated way. So while your "9" might now be classified correctly, the behaviour of the network on all the other images is likely to have completely changed in some hard-to-control way. That makes it difficult to see how to gradually modify the weights and biases so that the network gets closer to the desired behaviour. Perhaps there's some clever way of getting around this problem. But it's not immediately obvious how we can get a network of perceptrons to learn.

We can overcome this problem by introducing a new type of artificial neuron called a sigmoid neuron. Sigmoid neurons are similar to perceptrons, but modified so that small changes in their weights and bias cause only a small change in their output. That's the crucial fact which will allow a network of sigmoid neurons to learn.

The networks described are called feedforward neural networks. This means there are no loops in the network - information is always fed forward, never fed back. There are also other models of artificial neural networks in which feedback loops are possible. These models are called recurrent neural networks. The idea in these models is to have neurons which fire for some limited duration of time, before becoming quiescent. That firing can stimulate other neurons, which may fire a little while later, also for a limited duration. That causes still more neurons to fire, and so over time we get a cascade of neurons firing. Loops don't cause problems in such a model, since a neuron's output only affects its input at some later time, not instantaneously.

Deep Neural networks involve a network of layers which breaks down a very complicated question - does this image show a face or not - into very simple questions answerable at the level of single pixels. It does this through a series of many layers, with early layers answering very simple and specific questions about the input image, and later layers building up a hierarchy of ever more complex and abstract concepts. Networks with this kind of many-layer structure - two or more hidden layers - are called deep neural networks.

Since 2006, a set of techniques has been developed that enable learning in deep neural nets. These deep learning techniques are based on stochastic gradient descent and backpropagation, but also introduce new ideas. These techniques have enabled much deeper (and larger) networks to be trained - people now routinely train networks with 5 to 10 hidden layers. And, it turns out that these perform far better on many problems than shallow neural networks, i.e., networks with just a single hidden layer. The reason, of course, is the ability of deep nets to build up a complex hierarchy of concepts. It's a bit like the way conventional programming languages use modular design and ideas about abstraction to enable the creation of complex computer programs.


Denoise Techniques When Clean Signal and Pure Noise Are Available I have the clean version of the signal.
I can obtain the environmental noise.
I want to apply an effective denoising technique on a noisy signal (i.e., clean plus environmental noise).

Some observations:

The noise to signal ratio is extremely low.
The noise is spread across all frequencies in the frequency domain.

Are there any techniques for possibly taking advantage the noise and the clean signal towards denoising the noisy signal?

If the data is stationary and the noise is white then you should use the Wiener Filter.

If data isn't stationary you should look into the family of adaptive filters (LMS Filter and RLS Filter for starter).


downsampling an fft signal I am currently working with a set of wavelet transforms that I came across in a paper. The paper suggests that I first apply the FFT on the source image and then perform the downsampling. So it says something like:

Consider a source image of size 128*128.

P = fft(source);
P(1:64,1:64) = P(1:64,1:64) + P(65:128,65:128);
P(65:128,1:64) = P(65:128,1:64) + P(1:64,65:128);
P = P(1:128,1:64);

Well this solution definitely works and I get the results, but could anyone give a possible explanation as to why I am supposed to add in this manner?

Downsampling and spectral stretching is earlier to show in 1D. I'll also analyze with the Fourier transform and discrete-time Fourier transform, not a DFT. The results are the same except for aliasing in time.

To downsample by an integral multiple is to increase the sample period of a continuous-time signal $x(t)$ by an integer $M : T → M T$ . The replication period in the frequency domain is reduced by the same multiple. Rather than the first replica being centered at $\omega = 2\pi/T$ , the first replica of the down-sampled spectrum is centered at $\omega = 2\pi/MT$. Moreover, the amplitude is scaled by $1/M$. The figure illustrates $M = 2$. Sequencing the original series converts $x(nT ) → x[n]$. To downsample the sequence is to introduce a stride $M$ when reading off the values:

$$
x^{\downarrow} [n] = x[nM ]
$$

The effect on the angular sample frequency can be shown by a units analysis,

$$
\left(\frac{φ^↓\; \text{radians}}{1\; \text{sample}}\right)
=
\left(\frac{ω\; \text{radians}}{1\; \text{sec}}\right)
\left(\frac{T\; \text{sec}}{1\; \text{step}}\right)
\left(\frac{M\; \text{step}}{1\; \text{sample}}\right)
$$

or $φ^{\downarrow} = \omega M T$ . The range $ω ∈ [0, 2π/M T ]$ corresponds to $φ^↓ ∈ [0, 2π]$, thus $φ^↓$ is a stretched version of $ω$. The figure below compares $x[n]$ to $x^↓ [n]$ and $φ$ to $φ^↓$ . $x^↓ [n]$ has fewer samples than $x[n]$, and has a broader spectrum when plotted in $φ$ because the radians per sample has increased by $M$.

So in the 1D analogy to your case, refer to the continuous-frequency spectra $\omega$ above. Downsampling by 2:1 creates a spectral replica centered at $2\pi/MT$, $M=2$. That is what your first lines of code do (in 2D). Then, the spectrum is stretched so that the periodicity of the spectrum is restored to $2\pi$. This is shown to the right.


How to get the values for C1 and C2 in Chan-Vese Model? This question will consider the widely accepted segmentation algorithm of Chan-Vese, given in the following publications:

  T.F. Chan, L.A. Vese, An active contour model without edges, Lecture
  Notes in Computer Science, vol. 1682, pp. 141–151, 1999.

  T.F. Chan, L.A. Vese, Active Contours Without Edges, IEEE
  Transactions on Image Processing, vol. 10, no. 1, pp. 266–277, 2001

The main publication is accessible here:
http://www.math.ucla.edu/~lvese/PAPERS/IEEEIP2001.pdf

Equation 6 and 7 of this paper mentions $c_1$ and $c_2$ values. The energy function is also parameterized by these variables. My question involves the algorithm / method to compute these values. How can I do that?

The actual formulas are given on the publication, which I linked by editing your question. Besides, there are many publicly available implementations of Chan-Vese algorithm, including, but not limited to:

https://www.mathworks.com/matlabcentral/fileexchange/23445-chan-vese-active-contours-without-edges

https://github.com/kevin-keraudren/chanvese

http://www.ipol.im/pub/art/2012/g-cv/

https://www.mathworks.com/matlabcentral/fileexchange/19567-active-contour-segmentation

http://math.arizona.edu/~rcrandall/code.html

I think you could find all the details in those generously shared source codes.


Signal Spectrum Shaping I am looking for methods/algorithms for signal spectrum shaping. Let's say we have a white Gaussian noise and want to have a desired power-spectral density. The first solution that comes to mind is filtering. I was wondering if there are any other techinques.
I appreciate any hints.

but... "shaping the spectrum" is exactly the definition of a filter. I think you can justifiably call everything that shapes a spectrum a filter, so no, nothing but filtering shapes a spectrum.
Note that just because that's the most commonly analyzed case, filters are more than just linear filters.


Interpreting the cross correlation results I am new to this forum. I am just trying to understand how to calculate delay (in terms of samples) between two signals using cross-correlation function. I wrote a small Matlab code in order to understand this. But I am finding difficulty in understanding the results.
 clc;

 fm = 1000; //message signal
 fs = 10 * fm; //sampling frequency
 ts = 1 /fs;
 t = 0:ts:1-ts;

 x=sin(2*pi*fm*t);
 y=cos(2*pi*fm*t);

 [r,s] = xcorr(x,y);

 subplot(2,1,1);
 plot(r);
 subplot(2,1,2);
 plot(s);

The results obtained using this code are as follows.

From the figure, we can observe that the peak is occurring at around 10,000 sample. So, can 10,000 samples be considered as a delay (in samples) between two signals? And does the corresponding time delay is 10000 * 1/10000 = 1 sec in this case?

Yes and no. In principle, you can use the peak of your correlation function. However, it is not at 10000. The correlation function is symmetric around 0, so your peak is actually much closer to zero than you think. This is one of the reasons why xcorr returns two paramters, one for the lags at which the function is calculated. The correct way of plotting the the correlation function would thus be a plot(s,r).
Secondly, you should try zooming in a bit. Since your time signals have 10000 samples, your acf has ~20000, making it hard to see the exact maximum. If you do zoom in you'll notice it's not exactly at zero but two samples to the right.
Why two samples? Because you have 10 samples per period and the shift between sine and cosine is a quarter period, so it's actually 2.5 samples. You notice the problem, with the sampling you've chosen it is a little hard to measure the delay precisely.
Try reducing your fm to a small value, like 10, this makes the correlation function easier to interpret. Here is an example:

In this example, the maximum occurs at sample 247, which is equivalent to 247*ts = 24.7 milliseconds.
In case you wonder about the triangular shape, this has to do with the fact that while at lag 0 (in the middle) we can "compare" the entirety of the two signals whereas for larger lags, we need to truncate both signals to shorter versions. The larger the lag, the fewer samples we have to correlate. For this reason, we typically check the ACF closer to lag zero.


If group delay is a constant for all $\omega$, Does this system have a linear phase response? I want to know if a system has a linear phase response using frequency response of this system.

So I got phase response of it and also group delay.

and after that I got to know that if group delay is a constant for all ω, this system has a linear phase response.

so questions are

Group delay = $1/(1+\omega^2)$ . so $\omega=0$, group delay is $1$ and $\omega=1$ group delay is $1/2$. This is not a constant for all $\omega$ so I think this system has not a linear phase response. Is it correct?
when $H \left( e^{j\omega} \right)$ = $-e^{-j\omega} + 2 e^{-2j\omega} - e^{-3j\omega}$. Can I also know if this system has a linear phase response with same logic in question 1?

This is a homework type question, so I will only give hints and no solutions. You should know that the group delay is the negative derivative of the phase (with respect to frequency), so if the phase is a linear function of frequency, the group delay must be a constant. This should answer your first question.

Concerning your second question, why don't you just compute the phase of the given frequency response? It may be helpful to rewrite $H(e^{j\omega})$ as

$$H(e^{j\omega})=2e^{-2j\omega}\left(1-\frac{e^{j\omega}+e^{-j\omega}}{2}\right)\tag{1}$$

and to realize that the term in parentheses is real-valued.

EDIT:

The frequency response $H(e^{j\omega})$ can be written as

$$H(e^{j\omega})=M(\omega)e^{j\phi(\omega)}\tag{2}$$

where $M(\omega)$ is the magnitude, and $\phi(\omega)$ is the phase. Now compare $(1)$ and $(2)$ and figure out what $\phi(\omega)$ must be.


What's the difference between "Discrete Fourier Series" and "Discrete Fourier Transform"? I look at the equations for DFS (Discrete Fourier Series) and DFT (Discrete Fourier transform) and the only difference I notice is that one has a squiggle above the letter and the other doesn't. The equations for x[n] and X[k] are exactly the same except for the presence or absence of a squiggle.

Example:

DFS:  $\tilde{x}[n]$ / $\tilde{X}[k]$

compared to:

DTF:  $x[n]$ / $X[k]$

What's the difference between "Discrete Fourier Series" and "Discrete Fourier Transform"?

They are the same thing. Just the same.

But for some reasons, we want to call the first period of the DFS sequences $\tilde{x}[n]$ and $\tilde{X}[k]$, as the new DFT sequences $x[n]$, $X[k]$.

They are both periodic. The periodicity of DFS is natural extension of periodicity of complex exponential and periodicity of time sequence, but the periodicity of DFT is also forced by the modulus operator.

Be careful! DFT is also linked with finite length aperiodic sequences $x[n]$ and their DTFT via the frequency sampling relation; i.e., $X[k] = X(e^{j \frac{2\pi}{N}k})$.

Because of this latter relation of the DFT $X[k]$ to the DTFT $X(e^{j \omega})$ of aperiodic (and finite length) sequences $x[n]$, in order to remember and imply the inherent periodicity in DFT (due to its definition through DFS) we shall enforce the modulus operator into the arguments of DFT sequences.


sources for learning Image processing with MATLAB? I have Digital Image processing (DIP)  in my academics and I am looking for some reliable and daily updating sources for learning DIP Filters for Enhancements,color conversions, Filtering purposes through MATLAB , please suggest me some websites or tutorials which can help me. Thank you.

It does not update daily, but Steve Eddins' blog on image processing is an excellent source.

Steve is a co-author of Digital Image Processing using MATLAB, which is also a great source.

Another excellent source on image processing in general (not MATLAB-specific) are the lecture notes of Prof. Alan Peters from Vanderbuilt University.


Why is this FFT code not working to extract frequencies from a wav audio? I'm quite new to audio processing and i'm trying to extract frequencies from a wav file. As research papers stated i put the wav audio into a byte array and segmented and trying to apply FFT to get the frequencies from the audio. However the problem is even when i enter a audio track with only one note played (With no noise) i get a lot of frequencies as a result where as i'm supposed to get one.

My guess is something is wrong with the FFT code because the audio gets segmented properly.

Thanks everyone in advance,

My fft code

public void getFFT(double[] x, int a) {

double amp[] = new double[512];

double firstSum;
double secSum;
double twoPInjk;
double[][] s = new double[512][2];
for (int j = 0; j < 512; j++) {
    firstSum = 0;
    secSum = 0;
    FreqArr[j] = (j * sampleRate) / (512);
    for (int k = 0; k < 512; k++) {
        twoPInjk = ((2 * Math.PI) / 512) * (j * k);
        firstSum += x[k] * Math.cos(twoPInjk);
        secSum += x[k] * Math.sin(twoPInjk);
    }
    amp[j] = Math.abs(Math.sqrt(Math.pow(firstSum, 2) + Math.pow(secSum, 2)));
}
highest(amp, a);
}

public void highest(double[] amp, int a) {
String note;
double all = 0;
double max_magnitude = -1;
int max_index = -1;

for (int i = 0; i < (512); i++) {

    all = all + amp[i];

    if (amp[i] > max_magnitude) {
        max_magnitude = amp[i];
        max_index = i;
    }
}
double freq = 0;
int cnt = 0;

freq = max_index * (sampleRate / 2) / (512);
System.out.println("Maximum_index: " + max_index);
System.out.println("Sample rate: " + sampleRate);

  System.out.println(freq + " is the frequnecy");
}

my result (for the 1 second audio with only one note played)

Sample rate 44100 My data size 70376 My myBitsPerSample 16 My channels 2 My myByteRate 176400 My Length 0.7983219954648526 0.7983219954648526 44100.0 Max_index: 56 Sample rate: 44100.0 2411.71875 is the frequnecy Max_index: 256 Sample rate: 44100.0 11025.0 is the frequnecy Max_index: 256 Sample rate: 44100.0 11025.0 is the frequnecy Max_index: 129 Sample rate: 44100.0 5555.56640625 is the frequnecy Max_index: 511 Sample rate: 44100.0 22006.93359375 is the frequnecy Max_index: 256 Sample rate: 44100.0 11025.0 is the frequnecy Max_index: 152 Sample rate: 44100.0 6546.09375 is the frequnecy Max_index: 422 Sample rate: 44100.0 18174.0234375 is the frequnecy Max_index: 26 Sample rate: 44100.0 1119.7265625 is the frequnecy Max_index: 462 Sample rate: 44100.0

and it goes on.

Your DFT is working as expected at extracting frequency information.  One musical note at one pitch usually contains a lot of spectral frequencies: from harmonics and overtones plus broadband noise from any transients, such as note attack and decay.  Furthermore, the max magnitude frequency from an actual sound recording is often not the same as the heard note pitch frequency (one possibility is the magnitude peak is instead a strong overtone).

Instead, search for pitch detection/estimation algorithms.  The good (working, more reliable) ones will do a lot more than just an FFT peak magnitude search.


Autocorrelation for periodic signals Autocorrelation for power signals is defined by $$R_x(\tau)=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^Tx(t)x^*(t-\tau)dt\tag{1}$$ Is it true that for periodic signals $(1)$ can be computed by $$R_x(\tau)=\frac{1}{T_0}\int_{T_0}x(t)x^*(t-\tau)dt\tag{2}$$What are the other ways for computing $R_x(\tau)$? In this question, it's showed that $$|X(f)|^2 =\sum_{n=-\infty}^\infty |c_n|^2 \delta (f-nf_0)\tag{3}$$
But I don't think that's a valid derivation since, as it's pointed by Matt L, it involves product of  $\delta(f)$'s which is meaningless. So is it possible to prove $(3)$ by other means or it's a false statement? And, in the end, can we conclude that for a periodic power signal $$(\mathcal{F}R_x)(f) = S_x(f) = |X(f)|^2 = \sum_{n=-\infty}^\infty |c_n|^2 \delta (f-nf_0)$$
where $S_x(f)$ is PSD, holds?

The OP states that for a (deterministic) power signal $x(t)$, the autocorrelation function is defined as
$$R_x(\tau) = \lim_{T\to \infty}\frac{1}{2T}\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt\tag{1}$$ and then wonders whether in the case when $x(t)$ happens to be a periodic signal with period $T_0$, then it is true that
$$R_x(\tau) = \frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt ??\tag{2}$$
The answer is Yes, for the reasons given below.
Let us consider the integral $\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt$ in $(1)$ when $x(t)$ is periodic with period $T_0$ and $T = nT_0+\alpha$ where $0\leq \alpha < T_0$. Note that the integrand also has period $T_0$.

When $\alpha=0$, the range of integration is over $2n$ periods of the integrand and so
\begin{align}
\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt &= \int_{-nT_0}^{nT_0} x(t)x^*(t+\tau) \,\mathrm dt\\
&= 2n\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt.
\end{align}
It follows that
$$\frac{1}{2nT_0}\int_{-nT_0}^{nT_0} x(t)x^*(t+\tau) \,\mathrm dt = \frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt.$$
The astute reader will have seen that if $\alpha=\frac{T_0}{2}$, the range of integration is over an interval of length $(2n+1)T_0$ (that is, over $2n+1$ periods of the integrand) and so we have
$$\frac{1}{2(n+1)T_0}\int_{-(n+\frac 12)T_0}^{(n+\frac 12)T_0} x(t)x^*(t+\tau) \,\mathrm dt = \frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt$$
also.
For $\alpha \in \left(0,\frac 12\right)$ (or $\alpha \in \left(\frac 12,1\right)$, the integral is over an integer number of periods (say $N$) of the integrand plus an interval $\beta T_0$ of length less than $\frac 12T_0$ and so
\begin{align}
\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt &= N\int_0^{T_0} x(t)x^*(t+\tau)\,\mathrm dt + \Delta
\end{align}
where $\Delta$ denotes whatever the contribution is from the integration of $x(t)x^*(t+\tau)$ over that little interval of length $\beta T_0 < \frac 12 T_0$.
But then we have that
\begin{align}
\frac{1}{2T}\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt &= \frac{1}{(N+\beta)T_0}\left[N\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt + \Delta \right]\\
&= \frac{N}{N+\beta}\cdot \left[\frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt\right] + \frac{\Delta}{(N+\beta)T_0}.
\end{align}
But $(1)$ says that we need to take the limit of the above as $T \to \infty$ which we can do by letting $N\to\infty$ and so we see that in all cases,

If $x(t)$ is a deterministic periodic power signal, then its autocorrelation function
$$R_x(\tau) = \lim_{T\to \infty}\frac{1}{2T}\int_{-T}^T x(t)x^*(t+\tau) \,\mathrm dt\tag{1}$$
can be computed via
$$R_x(\tau) = \frac{1}{T_0}\int_0^{T_0} x(t)x^*(t+\tau) \,\mathrm dt, \tag{2}$$ which is what the OP asks for.

For the important special case when $x(t)$ is a real-valued function, Eq. $(2)$ simplifies to
$$R_x(\tau) = \frac{1}{T_0}\int_0^{T_0} x(t)x(t+\tau) \,\mathrm dt, \tag{3}$$
and this is the only case we consider below.

Note that so far, there has been no whiff of Fourier transforms or Fourier series or Dirac deltas or anything to do with the frequency domain, and so let us bring that stuff in.  Note that a real-valued deterministic periodic power signal $x(t)$ can be represented by a Fourier series
$$x(t) = \sum_{n=-\infty}^\infty c_n\exp(j2\pi nf_0t)\tag{4}$$ where $f_0 = T_0^{-1}$ and $c_{-n}=c_n^*$, and does not enjoy a Fourier transform in the classical sense (where Dirac deltas a.k.a. impulses are not allowed). However, $x(t)$ does have a Fourier transform $X(f)$ in the generalized sense where
$$X(f) = \sum_{n=-\infty}^\infty c_n\delta(f-nf_0).\tag{5}$$
In the classical theory, if $X(f)$ is the Fourier transform of $x(t)$, then $|X(f)|^2$ is the Fourier transform of its autocorrelation function $R_x(\tau)$. Does this result hold in the generalized sense too? That is, for a periodic signal $x(t)$ with Fourier series as in $(3)$ and periodic autocorrelation function $R_x(\tau)$ as in $(2)$, what is the Fourier series for $R_x(\tau)$ and the corresponding generalized Fourier transform of $R_x(\tau)$???
We have that $R_x(\tau)$ has Fourier series $\displaystyle \sum_{n=-\infty}^\infty a_n\exp(j2\pi nf_0\tau)$ where
\begin{align}a_n &= \frac{1}{T_0}\int_{0}^{T_0} R_x(\tau)\exp(-j2\pi nf_0\tau)\,\mathrm d\tau\\
&= \frac{1}{T_0}\int_{0}^{T_0} \left[\frac{1}{T_0}\int_{0}^{T_0} x(t)x(t+\tau)\,\mathrm dt\right] \exp(-j2\pi nf_0\tau)\,\mathrm d\tau\\
&= \frac{1}{T_0}\int_{0}^{T_0} x(t)\left[\frac{1}{T_0}\int_{0}^{T_0} x(t+\tau)\exp(-j2\pi nf_0\tau)\,\mathrm d\tau\right] \,\mathrm dt\\
&=\frac{1}{T_0}\int_{0}^{T_0} x(t)\left[\frac{1}{T_0}\int_{0}^{T_0} x(t+\tau)\exp(-j2\pi nf_0(t+\tau))\,\mathrm d\tau\right] \exp(j2\pi nf_0t)\,\mathrm dt\\
&=c_{n}\cdot \frac{1}{T_0}\int_{0}^{T_0} x(t) \exp(j2\pi nf_0t)\,\mathrm dt\\
&= c_n\cdot c_{-n}\\
&= |c_n|^2.
\end{align}
Thus, for a real-valued deterministic periodic power signal $x(t)$,

$$\text{If } x(t) = \sum_{n=-\infty}^\infty c_n\exp(j2\pi nf_0t), ~~\text{then }
R_x(\tau) = \sum_{n=-\infty}^\infty |c_n|^2\exp(j2\pi nf_0\tau). \tag{6}$$

Finally, what about Fourier transforms in the generalized sense?
Well, we have that
\begin{align}
\mathscr F[x(t)] &= \sum_{n=-\infty}^\infty c_n\delta(f-nf_0)\\
\mathscr F[R_x(\tau)] &= \sum_{n=-\infty}^\infty |c_n|^2\delta(f-nf_0)
\end{align}
and so it is not the case that the formula
$$\mathscr F[x(t)] = X(f) \implies \mathscr F[R_x(\tau)] = |X(f)|^2
\tag{7}$$
of classical Fourier analysis holds unless we are willing to assume that when we multiply out $\left|\sum_{n=-\infty}^\infty c_n\delta(f-nf_0)\right|^2$ into a doubly infinite sum over $m$ and $n$,
$$c_n \delta(f- nf_0)c_m^* \delta^*(f- mf_0) = \begin{cases}
|c_n|^2 \delta(f- nf_0), & \text{if}~m=n,\\
0, &\text{if}~m\neq n,
\end{cases}
\tag{8}$$
holds. In particular, we must treat $|c_n|^2\delta(f-nf_0)\delta^*(f-nf_0)$ as being equal to $|c_n|^2\delta(f-nf_0)$ as stated in my previous answer that the OP has denigrated as invalid and asserted that the claim that, in this context, the square of a Dirac delta should be treated as a Dirac delta, is meaningless because $\delta^2(x)$ is meaningless. The alternative is to insist that $(7)$ holds in classical Fourier analysis only and cannot be applied in the generalized Fourier analysis that allows for Dirac deltas. You pays your money and you takes your choice.....


FIR filter digital differentiator with low cutoff I'm trying to design a digital differentiator FIR Filter. It features a lowpass, such that above the cutoff frequency the amplification is very low. I get the coefficients by a linear program minimizing the chebychef error of desired and actual frequency response.

It works really well, but I cannot place the cutoff frequency below some 0.1*pi rad/sample. Small cutoff frequencies still have very steep rising amplitude responses in low frequencies and thus need a broad transition band.

The picture shows such a design and the very broad transition band. The red is the desired and blue obtained frequency response. I've weighted the bands accordingly. Also I'm not talking about bandpass, nor lowpass, I design a differentiator - thus the linear slew rate in low frequencies.

There are limits to the possible lowpass frequency, correct? How can I make the cutoff even smaller, or even better: why is this degradation happening?

I know, that the frequency response in my formulation has the form
$$
H(e^{j\omega}) = 2\sum_{k=0}^M j \,h(k)\, \sin(k \omega)
$$
where $M$ is $(N-1)/2$ with order $N$ filter. And thus the shape can be better traced by having longer filters. But the gain actually is very small.

Also I read, that with a derived then sampled Blackman Window (without control over cutoff frequency) one obtains a cutoff of around $\omega_C \approx 0.005$, while I struggle with $0.1$... I want to know why exactly.

This document suggests a method for a first order differentiator, where "only" the derivatives at $\omega = 0$ are matched to the ideal one. It results in an earlier drop off. However, as I understand it, this cannot be achieved with higher order differentiators, since second order is a quadratic function and I am not sure if a (basically taylor) approx in derivatives is sufficient for that. Let alone even higher orders.

What you want is what I call a Jekyl and Hyde design. The low pass is parsimoniously represented in terms of cosines and the derivative in sines. So let Jekyl be Jekyl and Hyde be Hyde.

The way to go very low would be a multirate approach.  Low pass and decimate for a few cascades until the high pass part of your differentiator is a reasonable fraction of the band.  If you need the original sample rate, upsample appropriately.  Your filter plot implies that you really aren't interested in most of your original band, so why retain it?


A system that perfoms Fourier Transform operation - is it an LTI system? If a system takes input as the time domain signal and outputs the frequency domain signal, is such a system an LTI system?
 For if the input time domain signal can be represented as a linear combination of sinusoids, then since the operation is transformation to the frequency domain, the outputs will also be the same Linear Combinations of the individual outputs (frequency domain representation of the individual sinusoids).. right? so we can say it's Linear?
 Since the system performs the same operation independent of time, it would be Time-Invariant too right?

The Fourier transform operator $\mathscr{F}$ is a linear one; i.e.,
$$\mathscr{F}\{x(t)\}=X(f) ~,~ \mathscr{F}\{y(t)\}=Y(f)  \implies \mathscr{F}\{\alpha x(t) + \beta y(t) \} = \alpha X(f) + \beta Y(f)$$ And therefore the system that implements it will be linear.

However time invariance, which is tested by the method

$$\mathcal{T}\{x(t)\}=y(t) \implies \mathcal{T}\{x(t-d)\}=y(t-d)$$

does not apply to the Fourier Transform operator, because the input and output functions do not belong to the same domain (namely time here)

However, the physical system which produces Fourier transform of the given signal $x(t)$ can be associated with a time-dependent Fourier transform which computes the function

$$X(\omega,t)= \mathcal{TF} \{ x(t) \} = \int_{-\infty}^{\infty} x(t+\tau) w(\tau) e^{-j \omega \tau} d\tau $$

Then it can be shown that this operator is time invariant in the sense that,

$$ \mathcal{TF} \{ x(t) \} = X(\omega,t) \implies \mathcal{TF} \{ x(t-d) \} = X(\omega,t-d)$$

Lets elaborate on this. Now lets define STFT as
$$X(\omega,t) = STFT \{ x(t) \} = \int_{-\infty}^{\infty} x(t+\tau) w(\tau) e^{- \omega \tau} d\tau $$ It's easy to show that $$X_2(\omega,t) = STFT \{ x(t-d) \} = \int_{-\infty}^{\infty} x(t-d+\tau) w(\tau) e^{- \omega \tau} d\tau = X(\omega,t-d) $$ Hence it's a shif-invariant operator.

EDIT: Below shown is a trivial redefinition of the above STFT, which does not confirm to the standard notation, hence may not apply. I don't delete it for now however. Shift-invariance is shown for the above definition only.

Now a similar form discussed in the provided link is obtained by a change of variable $t+\tau = \tau'$ on the STFT definition integral:

$$ X(\omega,t) =  \int_{-\infty}^{\infty} x(\tau') w(\tau'-t) e^{-j \omega (\tau'-t)} d \tau' $$
and replacing $\tau'$ with $\tau$ yields the equivalent form of STFT as
$$X(\omega,t) = STFT \{x(t)\}= e^{j \omega t} \int_{-\infty}^{\infty} x(\tau) w(\tau-t) e^{-j \omega \tau'} d \tau $$

Now lets show its time invariance by the input $x_2(t) = x(t-d)$ :
$$ X_2(\omega,t) = STFT \{ x_2(t) \} = STFT \{ x(t-d) \} = e^{j \omega t} \int x(\tau-d) w(\tau-t) e^{-j \omega \tau} d\tau $$ Make a change of variabes inside the integral operator as $\tau-d = \tau'$ yields

$$ X_2(\omega,t) = STFT \{ x(t-d) \} = e^{j \omega t} \int x(\tau') w(\tau'-(t-d)) e^{-j \omega (\tau'+d)} d\tau' $$ which simplifies by replacing $\tau'$ with $\tau$ and extracting $e^{-j\omega d}$ outside of the integral

$$ X_2(\omega,t) = STFT \{ x(t-d) \} = e^{j \omega t} e^{-j \omega d} \int x(\tau) w(\tau-(t-d)) e^{-j \omega \tau} d\tau $$

$$ X_2(\omega,t) = STFT \{ x(t-d) \} = e^{j \omega (t-d)}  \int x(\tau) w(\tau-(t-d)) e^{-j \omega \tau} d\tau = X(\omega,t-d)$$

Which proves that this form is also shift-invariant on the $t$ variable.


Periodicity of Fourier series's coefficients Why exactly continious Fourier-series's coefficients aren't periodic like coefficients of discrete Fourier series (DFS)? $e^{-j2\pi}$ is periodic in both sequences.

The reason why Fourier series coefficients of continuous functions are generally not periodic is because of the continuity of the function. The discrete Fourier series coefficients are periodic because the analyzed signals are discrete. Note the duality relationship of the Fourier transform

$$\text{periodicity}\Longleftrightarrow\text{discreteness} $$

Periodicity in one domain leads to discreteness in the other domain. E.g., a periodic function has only discrete frequency components (as shown by its Fourier series). If these discrete frequency components are periodic, then, consequently, the periodic function must be discrete, i.e., it must be non-zero only at discrete times.

In his answer, Fat32 gave an example of a periodic signal with periodic Fourier series coefficients. The most general form of such a signal is

$$x(t)=\sum_{m=-\infty}^{\infty}b_m\delta\left(t-m\frac{T}{N}\right)\tag{1}$$

where the coefficients $b_m$ are essentially the discrete Fourier series coefficients of the $N$-periodic Fourier coefficients $a_k$ of the $T$-periodic signal $x(t)$:

$$b_m=\frac{T}{N}\sum_{k=0}^{N-1}a_ke^{j2\pi km/N}\tag{2}$$
$$x(t)=\sum_{k=-\infty}^{\infty}a_ke^{j2\pi kt/T}\tag{3}$$

As mentioned above, the signal $x(t)$ is of course discrete in the sense that it is non-zero only at discrete time instances $t_m=mT/N$. Any periodic signal $x(t)$ with $N$-periodic Fourier coefficients $a_k$ must be of the form $(1)$.

A derivation of Equations $(1)$ and $(2)$ can be found in this answer.


Beginner's book in signal processing with practical examples on fault detection in electrical motors I'm an electrical engineering graduate and I'm working on condition monitoring for electrical motors and generators. My goal is to improve reliability of electric rotary machines in such way, that any ongoing faults are discovered before they progress in major failures. Main role in this domain has signal processing (and also statistical analysis along with machine learning), which I would use to develop condition monitoring algorithms. Unfortunately signal processing wasn't included in my education program, so I have to learn this on my own.

So what I think that I need, are some sources which outline  basic concepts, typical tools (FFT, STFT, wavelets, filters...) and problem-solving examples (noise removal, hidden pattern extraction, ripple reduction,...)  with only so much math involved that the main idea can be grasped, and that implementation is possible, which would be demonstrated in lots of MATLAB examples in the context of electrical rotary machines. I truly don't want a source for signal processing, which is aimed for mathematicians or information engineers, but rather sth applicable for electrical engineers.

When I had a look at rotating machinery, the best reference I could find is Bob Randall's Frequency Analysis.  This was generated in conjunction with Brüel & Kjær as they sold lots of nice (and expensive) vibration measuring / monitoring equipment.  A short Google search shows that B&K have made the PDF of this book available for free download from their website.

The book is targeted at mechanical engineers, who (back in the day) were never exposed to the Fourier transform.


Matlab: Detect Peaks and Signal Jump I want to detect 6 peaks in a signal. I am open to any suggestions. There is one problem where my signal sometimes jumps high and stays higher than the rest of the signal. I don't know why this happens but would like a solution to lower the signal back to the points where the rest of the signal was lower, which might help detect all 6 peaks. I attached a picture.

I think the easiest option is to detect the massive sudden increase/decrease.
Try calculating the difference between samples. If this difference is too big (positive or negative), we substract this difference on all future measured samples.

Create a variable for this offset. start with an offset of 0, when de absolute value of the difference between samples is too big, add this (non absolute) difference to the offset, then all next samples are $x_{new}=x-offset$.


What is the bandwidth of human speech? What is the actual frecuency bandwidth of human speech.
I know there is variation for gender, age and so, but what is the general range... I seem to find different ranges.

Like here http://recognize-speech.com/preprocessing

Is the range defined between 100 hz 8khz

Telephones are made to work in the ranges 300 Hz to 3400 Hz, according to Wiki: https://en.wikipedia.org/wiki/Voice_frequency

The hearing is in range

That really depends on your requirements and on your exact definition of bandwidth.

Telephone range (300 Hz-3.4 kHz) has okay intelligibility but the quality of the voice is fairly compromised. For example it's hard to tell the difference between "s" and "f" sounds.

80 Hz to 8 kHz gives you very good quality, although it's still possible for the human voice to make sounds outside of this range. Someone like James Earl Jones can easily go don to 50 Hz and specially trained singers can go a lot lower. Some fun examples at https://www.youtube.com/watch?v=lkdMl1ITw5M

If you want all of it, you should go for the range of hearing which is about 20 Hz to 20 kHz.


Convolution including $\delta(t-5)$ I know the two properties of convolution that are related to my question

$\quad x(t)*\delta(t)=x(t) $
$\quad x(t)*\delta(t-t_0)=x(t-t_0)$

But my question is, how can I use those two to calculate $$ y(t)=x(5t)*\delta(t-5) $$

HINT:

Set $f(t)=x(5t)$ and use your rule number 2: $$f(t)\star \delta(t-t_0)=f(t-t_0)=\ldots$$


FFT Data Array Indexing I've been Googling "FIR filters" and "FFT convolution" etc. and I ran across something interesting.

Array Index vs. Pointers:
http://www.bores.com/courses/intro/program/7_array.htm
http://www.dspguru.com/dsp/tricks/fir-filtering-in-c

So, here's my question:
Why is almost every "Fast!!" Fourier Transform function on the web using Array Indexing?
(I just thought I would throw it out there.)

First link is broken, the second link is about unrolling loops, so it's not entirely clear (to me) what exactly you mean by "array indexing". Anyway most of this heavily depends on the instruction set of the processor. For example and FIR filter would have to do something like this (in the innner loop)

read filter coefficient
read delayed input sample (state variable)
multiply the two
add to accumulator
decrement coefficient pointer (or decrement coefficient read index)
increment state pointer (or increment state read index)
warp state pointer if needed (these tend to be circular buffers)
decrement loop counter, test for 0
jump to beginning of loop or past the loop depending on previous test

So that's a lot of stuff to do. Especially the "wrapping" part can be quite expensive as it includes a test for wrap around and a conditional which can break the pipeline on many processors.

Anyway on some processors, like a risc, or i86 this can easily take 10-20 cycles, on a dedicated DSP like an Analog devices Sharc or a TI TMS320 these take about half a cycle (can do two in a single CPU clock), so there is a huge dependency on instruction set efficiency. The tricks the second long talks about can help with avoiding the unwarpping part and also reduces the overhead for looping since you only do it once every 4 coefficients.


Building a Signal Model for Samples from a Sensor I have a signal I'm getting from a sensor (shown below in the photo of an excel plot) and I need to process it to get it to look something like the red line I've overdrawn on it.  Sort of a moving average type process that doesn't fall too quickly unless the overall signal is actually falling.

Anyone have any ideas as to what the best method(algo) would be to approach this problem?

Any suggestions on directions I could look towards would be helpful as well, it's been years since I took signal processing in school and I don't remember what methods are even out there that could be tried.

Thanks!

Jonathan

A simple exponential average has the following difference equation:

y(n) = (1-alpha) * x(n) + alpha * y(n-1)

To get a different rise and fall rate, you can dynamically choose alpha depending on whether x(n) is larger than y(n-1). For example, you could build have alpha = 0 when x(n) > y(n-1), and alpha = 0.99 otherwise. This could accomplish the "peak hold" goal, and the second value of alpha would let you control the decay after each peak.