Freeform band-limited waveforms I'm ready to dip my toes in the scary waters of non-aliasing waveforms and wavetables. Ideally, what I'd like to do is create a (naive) wavetable/waveform in memory and subsequently generate band-limited versions of it with a reusable procedure.

Is that possible, and what would be the best starting point for my research?

this paper was done long before i had MATLAB.  the drawings are poor, but the math (at least in this revision, which is what you should use) is spot on.

send me an email address (to my audioimagination.com, not the wavemechanics.com on the paper) and i will send you a very short C file that shows how to generate and crossfade the wavetables in the synthesis process.

about analyzing a sample of a note (like from a wav file), you need to have a good sliding pitch-detector and a good interpolator (like maybe a good Kaiser-windowed $\operatorname{sinc}(\cdot)$-based interpolator) to extract single-cycle waveforms out of the note and put them into wavetables.

or you can create wavetables in memory from the POV of additive synthesis, simply by adding up the sinusoids, but it might be easier with the inverse FFT (and this is when knowledge of the periodicity of the DFT can be used to advantage).  i think your life will be easiest that the active wavetables in memory that are being used in the real-time synthesizer, those wavetables should all be the same power of 2 (and a nice big power of 2, like 1024 or 2048 or 4096).  notes that are higher in pitch will have fewer harmonics than notes lower in pitch, but will have the same number of samples per wavetable.  (and if the number of samples greatly exceeds twice the number of harmonics, then you can get away with linear interpolation during synthesis.)

wavetables that are being stored and not active, need not have so many samples.  as long as the number of samples in the periodic wavetable is more than twice the highest harmonic number, then, at least theoretically, it is enough.  but when a wavetable is made "active" and loaded into memory, it should be expanded and interpolated out to be whatever the large power of 2 is.

oh, and Nigel Redmon has code somewhere.


Minimum bit-width for fixed point FFT outputs I was wondering if there is a formula to determine what the best case (minimum) bit-width should be for the imaginary and real outputs of an fixed-point FFT, assuming an input bit width of $M$ ($M$ bits for the real and $M$ bits for the imaginary portion of the input). Best case meaning there is no possibility of overflow, and I maintain the same number of fractional bits as the input.

I can sort of see that for a $N$-point FFT, there will be $\log_2(N)$ butterfly "passes" or "columns", and always the twiddles have complex magnitude of $1.0$, so a simple assumption would be: $M + \log_2(N)$. I think this is treating it as a worst-case adder assumption. But I know that the twiddles vary a lot after the first column, so this is probably not "best case" analysis.

it really depends on how the DFT is defined.  usually we define the DFT and inverse DFT as:

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \ e^{-j 2 \pi \tfrac{nk}{N}} $$

$$ x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] \ e^{+j 2 \pi \tfrac{nk}{N}} $$

but they could just as well be defined as

$$ X[k] = \frac{1}{N} \sum\limits_{n=0}^{N-1} x[n] \ e^{-j 2 \pi \tfrac{nk}{N}} $$

$$ x[n] = \sum\limits_{k=0}^{N-1} X[k] \ e^{+j 2 \pi \tfrac{nk}{N}} $$

all of the theorems work exactly the same, except there will be a scaling difference in the convolution theorem and Parseval's theorem.

and, my suggestion, for fixed-point arithmetic is splitting the difference:

$$ X[k] = \frac{1}{\sqrt{N}} \sum\limits_{n=0}^{N-1} x[n] \ e^{-j 2 \pi \tfrac{nk}{N}} $$

$$ x[n] = \frac{1}{\sqrt{N}} \sum\limits_{k=0}^{N-1} X[k] \ e^{+j 2 \pi \tfrac{nk}{N}} $$

now, let's assume for shits and grins that $N$ is a power of 2, so there will be $\log_2(N)$ passes in the fixed-point FFT code.

in the first definition (at top), in the inverse FFT, there will be a division by 2 (an arithmetic shift right) for the results of the FFT butterflies in each inverse FFT pass (and none of this in the forward FFT).  $X[k]$ will have a mean square that is $N$ times bigger than the mean square of $x[n]$.  so $x[n]$ might be expected to padded with $\log_2(N)$ extended sign bits on the left to keep $X[k]$ from getting too big for the word width.

in the second definition (in the middle), the roles are reversed and in the forward FFT, there will be a division by 2 for every FFT butterfly.  so, since the magnitude of $x[n]$ is expected to grow for each FFT pass, then $X[k]$ should be padded with $\log_2(N)$ guard bits on left.

in the third definition, $X[k]$ and $x[n]$ both have the same mean square or the same energy.  and in both the forward and inverse FFT, for every alternate FFT pass, the outputs of the butterflies are divided by 2.  now having the same mean square does not guarantee that word width limits are not exceeded.  if the words going into the first FFT (i'll assume the first FFT is the forward FFT) are padded with $\tfrac12 \log_2(N)$ guard bits, then i would not expect any overflow.  and if the processing in the frequency domain preserves the mean-square of $X[k]$, then i would not expect any overflow in the inverse FFT.

one totally different trick, if you have the code to efficiently do this, is to do "block-floating-point" arithmetic for the FFT.  but that might not be a very fast FFT because it requires testing the magnitude of the outputs of each butterfly.  the Mot DSP56002 and later chips (from the 1990s) had a sticky-bit in the status register that automatically detected this bit-width growth, and a selectable "divide-by-2" mode on the outputs, which both facilitated a very efficient block-floating-point FFT in the fixed-point processor.


Circular Convolution of Length L of Sequences of Length Greater Than L I'm trying to understand how may I obtain the circular convolution of length L when the sequences I'm trying to convolve are of length greater than L.

For instance this Matlab code using sequences of length 5:

c = cconv([0,0.5,1,1,0.5],[0,0.5,1,1,0.5],4);

yields a 4 point sequence

2.25  2.5  2.25  2.0

and this Matlab code

cconv([0,0.5,1,1,0.5],[0,0.5,1,1,0.5],3)

yields a 3 point sequence

3 3 3

I know that circular convolution may be seen as linear convolution with aliasing and that if I perform the linear convolution and build overlapping periods spaced every N samples I will get the N point circular convolution but is there a way to obtain circular convolution without having to compute first linear convolution? I'm thinking about the classical way of doing convolution, where I flip one of the sequences in time and shift it and for every shift I multiply the 2 sequences and add to obtain the output sample for that shift

Circular convolution is a very close relative to linear one. Both of them are computed based on the same principles. However circular convolution has a periodic result (based on the fact that it's inputs are periodic), hence its computation can be affected by this.

When implementing a circular convolution there is no loss of efficiency if you actually use a linear convolution underhood (vice versa) as their computational efficiencies are very close, except if you are dealing with architecture optimized code.

On the other hand you may be asking for the penalty paid for computing a longer than necessary linear convolution in order to implement an $L$ point circular convolution of $N$ point periodic sequences when $L < N$. In such a case you can improve the efficiency by the following chopping method (borrowed from @DilipSarwate 's a relevant answer ) which I give here a basic MATLAB implementation from which you can see that the linear convolution is performed at the shorter length $L$, instead of the longer one $N$.

Below gives you a way of computing $L$ point circular convolution of $N$ point sequences where $ L < N$.

clc; clear all; close all;

% S0 - Define the Signals
% Note this is meaningful only when L < signal lengths
L = 3;                 % circular convolution length that we want to compute
x = [0 0.5 1 1 0.5];   % signal length is greater than L
h = x;

xL = zeros(1,L);       % these are for short, chopped signals.
hL = zeros(1,L);
for k=1:L              % compute the CHOPPED signals...
    xL(k) = sum( x(k:L:end) );
    hL(k) = sum( h(k:L:end) );
end
yN = conv(xL,hL);     % compute the linear convolution at the chopped length
yL = yN(1:N) + [yN(L+1:end) 0] ; post-process the result.


Processing line scanner data I have a line scanner output that is scanning the geometry of a surface. Typically the output looks something like this:

Where the x and y axis refers to the coordinate system in mm. (ignore the red dots). This data is essentially a 2d point cloud.

So what i am trying to achieve is to detect the corner points of this scan. I try to start using image processing techniques, for example using a harris corner detector.

However, the issue with this is that this data has a resolution of up to 0.01um vs a possible area size of 140mm X 100mm of all possible data points. If i were to convert this to 2d image, it would be a very large file, with most of the pixels being empty.

My question is that is there a better way to go about processing this type of data; or is there an alternative way to apply techniques such as the harris corner without resorting to conversion into an image.

A first suggestion (without being sure to understand the question).

After sorting on X, try a line fit in a short sliding window (say 4 units long) and observe the correlation coefficient.

In smooth areas, the fit will be excellent. It will degrade at corners. Looking at the plot of the correlation, you should be able to set a suitable detection threshold.


LTI system phase response given $z$-transform I have been given this question

  \begin{equation}
H\left(z\right)\:=\:\frac{1}{6}\left(1+z^{-2}\right)^6
\end{equation}

  (a) Compute and plot the phase response of the system.

  (b) Determine analytically $\angle H(e^{j\omega})$ and use the formula obtained to compute and plot the phase response.

  (c) Compute and plot the phase response using the function freqz.

First of all how to input the expression into MATLAB?
I tried

Hz1 = [1 0 1];
Hz11 = Hz1.^6

but this does not work. Any clues?

Answer to (b):

$(1+z^{-2})$ is a Type I generalized linear phase filter with group delay 1 (i.e. the phase response is $\angle H(e^{j\theta})=-\theta$ with a jump of $\pi$ at $\theta = \pm \pi/2$).

6 of those filters in cascade will have group delay 6 (phase response $\angle H(e^{j\theta})=-6\theta$).


Viterbi Decoder - How does inverting the input bits affect the result? Does converting zeros into ones and ones into zeros in the input stream affect the decoding result? Empirically I have seen that it does not affect, decoding $[0 1 0 1 \ldots]$ and $[1 0 1 0 \ldots]$ produce the same output stream. But I am not sure if it is true or not, and I am looking for a theoretical proof.

Why I am doing this is because the convolutionally encoded stream (rate $\frac{1}{2}$ and constraint length 7) is received without the phase information, thus a conversion is needed. I am trying to determine in which step I should convert the data.

I would appreciate if you could share any knowledge.

Convolutional codes are linear in the Galois{2} (i.e. binary) field.  We know this because the convolution operation is linear.  The fact that Viterbi codes are linear has a few implications.  First, if you input all 0's you will get all 0's out.  Second, if you add two inputs together the output will be the sum of their individual outputs.

The last piece of information that we need for a very informal "proof" is that when you input all 1's the output is all 1's (EDIT: this is not always true- see the edit at the bottom of the answer).  We see that this must be so because the encoder state is constant (you are, after all, inputting all 1's) so the output must be constant.  The output cannot be 0's because that would match the all 0's case which would make the code non-linear (other than in the trivial case where the encoder always outputs 0 no matter what the input is).

Now that we've established that inputting all 1's causes the output to be all 1's we can prove that yes, inverting the input is equivalent to inverting the output.

$$
Input => Output\\
Input + Ones => Output + Ones
$$

If we assume the first statement then the second statement follows due to linearity and all 1's producing all 1's, with the understanding that binary inversion is equivalent to adding 1.  0 + 1 = 1, 1 + 1 = 2 = 0.

EDIT: Dilip is correct.  For some generator polynomials an all 1's input will produce an all 1's output and so inverting the input is equivalent to inverting the output.  For some generator polynomials, though, an all 1's input will not produce an all 1's output, so inverting the input is not equivalent to inverting the output.  The generator polynomials that do produce an all 1's output have odd Hamming weights for both polynomials.  The reason that Dilip's polynomial does not meet the criteria is that one of the polynomials ($1 + x^2$) has an even Hamming weight.

EDIT 2: Dilip makes a valid point that even polynomials with odd Hamming weights will produce some zeros at the beginning while the zeros in the state memory are shifted out and replaced with ones.  This is relatively brief, but it does mean that the inversion equivalence does not hold during that time.  It does hold during the "steady state" behavior.


Tips for improving pitch detection I'm working on a simple web app that allows the user to tune his/her guitar. I'm a real beginner in signal processing, so please don't judge me too harshly if my question is inappropriate.
So, I managed to get the fundamental frequency using an FFT algorithm and at this point the application is somehow functional. However, there is room for improvement. Right now, I send raw PCM to the FFT algorithm, but I was thinking that maybe there are some pre/post algorithms/filters that may improve the detection. Can you suggest any?
My main problem is that when it detects a certain frequency, it shows that frequency for 1-2sec and then jumps to other random frequencies and comes back again and so on, even if the sound is continuous.
I'm also interested in any other type of optimization if anyone has experience with such things.

I'm guessing the other frequencies it gets are harmonics of the fundamental?  Like you're playing 100 Hz and it picks out 200 Hz or 300 Hz instead?  First, you should limit your search space to the frequencies that a guitar is likely to be.  Find the highest fundamental you're likely to need and limit to that.

Autocorrelation will work better than FFT at finding the fundamental, if the fundamental is lower in amplitude than the harmonics (or missing altogether, but that's not an issue with guitar):

You can also try weighting the lower frequencies to emphasize the fundamental and minimize harmonics, or use a peak-picking algorithm like this and then just choose the lowest in frequency.

Also, you should be windowing your signal before applying the FFT.  You just multiply it by a window function, which tapers off the beginning and end of the waveform to make the frequency spectrum cleaner.  Then you get tall narrow spikes for frequency components instead of broad ones.

You can also use interpolation to get a more accurate peak.  Take the log of the spectrum, then fit a parabola to the peak and the two neighboring points, and find the parabola's true peak.  You might not need this much accuracy, though.

Here is my example Python code for all of this.


How to choose pilot symbols count in OFDM? I am from an Instrumentation background and recently have been working on an OFDM project. With my research, I was unable to find the solution to the below question.

How to choose the Pilot Carriers/pilot symbols? Suppose from one example which I am referring they have chosen pilot carriers as 8 counts for the  Subcarrier count of 64 and with Cyclic prefix is 16. and in the second example pilot carriers = 346 for the subcarrier count = 224 with CP = 16. I am really confused about how to choose the pilot count?

It is a tradeoff. The more pilot subcarriers, the better then channel estimation, but the fewer subcarriers usable for data. Simulations and testing of performance can help in making the choice.

You probably have some numbers wrong. If you have 224 subcarriers, you cannot have 346 of them as pilot carriers.


What is the physical interpretation of the absolute value of a fourier transformed signal, $\left| F(t)\right|$? If one has some oscillating voltage signal, for example:
$$V(t) = V_{max}\cos(2 \pi \nu_{0}t) e^{-\gamma t}$$
and you take the Fourier transform of this in the usual way to get:
$$\hat{V}(\nu) = V_{max}\int^{+\infty}_{0} \cos(2 \pi \nu_{0}t) e^{-\gamma t} \ dt$$
Now I often see in literature the absolute value squared taken of the result, so $\left|\hat{V}(\nu) \right|^{2}$ the physical interpretation of which it is the energy per unit frequency.

My question is, if $\left|\hat{V}(\nu) \right|^{2}$  is energy per unit frequency, then what is the interpretation of $\left|\hat{V}(\nu) \right|$?

This would be the so called magnitude spectrum, its unit is the same as the signal it's generated from, in the case of voltage its interpretation is "voltage  by time over frequency" as opposed to "power by time over frequency" of the squared magnitude spectrum.
In the ususal logarithmic representation, the two are proportional by a factor of 2 and essentially give you the same information.

The $\cos$ function in your example is a so called power signal, not an energy signal, as it is infinitly long it has infinite energy, so to assess its energy in a practical way, one refers to its power (energy per time).


Is the incoming sample in a standard DPCM already quantized? I began learning about DPCM and i've got one issue I am not sure about:

Is the incoming sample m[k] (see image for reference) already quantizied or not? I think it is but with an (assumed) very fine quanitization. So the Quantiser block is doing a requantization with a more coarse quantization, right?

No it is assumed that m[k] is not quantized, as the main purpose of DPCM is to quantize the difference between samples rather than quantize the whole samples.

ofcourse if m[k] is an image it must been quantized as the video sensor in any digital camera is a digital sensor by definition so it is doing some means of quantization.


Java - Convert Image to Binary (Two Colours/Black & White) Is it possible to convert a .jpg image to binary using Java?

I have asked a similar question on StackOverflow but thought that maybe this may be a good place to ask, too.

I am currently able to convert an image to negative in Java, but am unable as of yet to figure out a way of converting it into binary.

Any suggestions - and even sample code - would be amazing.

Many thanks.

Converting it to "negative" will give you single channel result (0-255). To convert this image to binary, you need thresholding. This threshold value will have a HIGH (0-255) and LOW (0-255). These values are dependent on your application (what you wish to segment). For instance, you could setup an empty image structure and set all values to '0'. Set '1' to channel values that fall between HIGH and LOW. However, remember that performing naive segmentation like this may not give you good results.


Initialize a Room Impulse Response using reverberation time(T60) I am doing Speech dereverberation using Non-negative matrix factorization. I need to initialize an impulse response($H$) such that its reverberation time ($t_{60}$) is $300$ ms. Let the length of $H$ be $10$.

I thought of doing that by identically initializing each row of $H$ using a linearly decaying envelope. Somewhat like this

x=1:10;
h = exp(-x/2);
H = repmat(h,600,1);

This will give me a $H$ of dimension $600$ X $10$. But, I don't understand how to use $t_{60}$ for the initialization in MATLAB. Please tell me.

A good method is to generate white noise apply an exponential envelope with the reverb time. In general you will need lots of samples. Typically 10s of thousands. If you want to get fancier you can window out the front part and put in a direct sound spike and a few discrete reflections. Below an example of how to create a simulated room impulse response with a given reverb time and sample rate

%% create an room impulse response with a given t60
fs = 44100; % sample rate in Hz
t60 = 0.35; % reverb time in seconds
% 1. Detemine the length of the impulse response. It's an infinite response
%    so some truncation is neccessary. A good starting point 1.5 times the
%    t60 which will result in 90 dB of dynamic range
n = round(fs*1.5*t60);
t = (0:n-1)'/fs;  % time vector
% 2. initialize to white noise, gaussian distributed
h = randn(n,1);
% 3. Calculate the decay. During the reverb time the envelope decays by 60
%    dB, so we have exp(decay*t60) = 1e-3; We can solve to
decay = log(1e-3)/t60;
% 4. Apply the envelope
h = h.*exp(decay.*t);
figure(1); clf
plot(t,h);
xlabel('Time in s');

%% we can calcuale the energy decay curve through reverse integration
decayCurveInDB = 10*log10(flipud(cumsum(flipud(h.^2))));
% normalize to 0 dB
decayCurveInDB = decayCurveInDB-decayCurveInDB(1);
% and plot it
figure(2);clf;
h = plot(t,decayCurveInDB);
set(h,'Linewidth',2);
xlabel('Time in s');
ylabel('Decay in dB');
% mark the t60 spot
hold on
set(gca,'ylim',[-100,0]);
y = get(gca,'ylim');
h = plot(t60*[1 1],y,'r');
set(h,'Linewidth',2);
grid('on');


Finding the Sections of Signals that are the Most Similar I have two discrete signals, $A$ and $B$ (both with $L$ values). I want to look for the section (of a certain length $R$) on each signal ($[A_1:A_1+R]$ and $[B_1:B_1+R]$) where they are the most similar (using some similarity metric, say Manhattan distance) (technically, I just need to know how similar the most similar sections of the two signals). Right now, I have a double loop that compares each section of signal $A$ to each section of signal $B$, and choose the $A_1$ and $B_1$ that lead to the most similar sections. However, that gets really slow as I scale up. Is there some other, more efficient way of doing this?

The best way to achieve this, to my knowledge, is a matched filter or cross-correlation. This will tell you not only how well $A$ and $B$ match but will also give you a time offset for it. Just make sure to normalized based on the energy of both signals so you aren't getting false positives if that's a major concern.
If you're looking for a more strict time match (the signal similarity needs to be in-time) then you only want to look at $R[0]$ where $R[n] = A[n] \star B[n]$.
Technically, cross-correlation of non Wide Sense Stationary signals is two-dimensional. So if you're interested in the particular time when the signals best match up with no offset, your process may look something like this over $Z$ trials.
$$
\begin{align}
R[n,m] &= E\{A[n]B^*[m]\}\\
m &= n\\
R[n] &= E\{A[n]B^*[n]\}\\
R[n] &= \frac{1}{Z}\sum_{\zeta=0}^{Z-1}{A_\zeta[n]B_\zeta[n]}
\end{align}
$$
If you clarify what sort of signals these are and what about their matching is important, I think it would help narrow down what solution would best fit.


How to convert sampled real audio signal to it's IQ components for IQ modulation? So my question is this:
How to convert sampled real audio signal to it's IQ components for IQ modulation?

For instance:  WHen I have my DSP or FPGA receive a audio signal through a ADC, which is just the RAW A/D value.  Now I want to run this through a DAC to IQ modulator to get let FM modulation.
How would I convert the raw A/D value into IQ?  How do i break this signal down to In-phase and quadrature component? Do I do FM modulation in the DSP/FPGA and get the IQ from that or just split the audio baseband into it's IQ components?

Added:

I'm in the process of converting FSK transmitter from older obsolete design to newer  design.  So as I'm studying IQ, i was wondering what would happen if I  wanted to send a analog sine wave through the IQ modulator, like LTC5599.  How would I convert a analog sine wave to I and Q send it as input to the LTC5599?

I’ll answer the question I think you’re asking here:

If you have a real valued signal that came right out of an A/D and you’d like to recover signal phase, you can use the Hilbert transform to generate a version of the signal that is in quadrature with your A/D output. In this case, your original real valued signal would be the I component (in-phase) and the output of the Hilbert transform multiplied by j, the imaginary number, would be your Q component (quadrature).

Now if you have some modulation you want to deal with, you would do the IQ recovery first like I just described; once you have that IQ signal you’re pretty much completely characterized the signal, so you can do just about whatever you’d like with it while preserving phase.


Need to calculate $\frac{E\left[|x(n)|^4\right]}{E\left[x^2(n)\right]}$ I need to calculate $\frac{E\left[|x(n)|^4\right]}{E\left[x^2(n)\right]}$ in real time, where $x(n)$ is complex.

For $E\left[|x(n)|^4\right]$ I plan to take the absolute, raise it to the power $4$ and do exponential averaging with weights $0.999$ and $0.001$. I have no clue on how to calculate $E\left[x^2(n)\right]$. I know that this gives the total signal power. Will this be equivalent to $E\left[|x(n)|^2\right]$?

Thanks, JK

The most efficient way that I can think of to do it would be as follows:

For each input sample $x[n]$, calculate $|x[n]|^2$. This is more efficient than calculating $|x[n]|$, because you don't have to do the square root:

$$
|x[n]| = \sqrt{\text{Re}\{x[n]\}^2+\text{Im}\{x[n]\}^2} \\
|x[n]|^2 = \text{Re}\{x[n]\}^2+\text{Im}\{x[n]\}^2
$$

Average the instantaneous power estimate $|x[n]|^2$ as needed for your application (you mentioned exponential averaging).
Note that you can also get $|x[n]|^4$ efficiently using this method; you just square the measurement of $|x[n]|^2$ that you made above, then average again as desired.


inverse Fixed point number (number^-1) Could someone please provide a method for inverting random fixed point numbers.
I have as an input both positive and negative floating point numbers that should be converted into fixed point and then be inverted.

An example could be 1/0.2453 where floating-point 0.2453 is converted into a Q18 fixed point format and then inverted to 4.07664.
Thanks a lot!!

the first thing that one must understand when doing fixed-point arithmetic is that, at the basic level, it is integer arithmetic with some scaling factors applied.  perhaps just one scaling factor applied.  this scaling factor is directly related to the position of the binary point.  unless the chip is a fixed-point DSP or similar, the binary point is implied, not explicit.  it, and the associated scaling factor, exists only in the mind of the programmer.  the CPU only thinks it is working on integers.

as an example an unsigned, fixed-point, 16-bit format, that i would call "1.15" or, perhaps "Q1.15" (i think the "Q" is superfluous and notation like "Q16" does not supply enough information) has 1 bit to the left of the binary point and 15 bits to the right.  but the 16 bits, sitting all by themselves, is an integer.  the range of that unsigned integer, $n$, is $0 \le n \le 2^{16}-1$.  but the range of the unsigned fixed-point number implied $\hat{n}$ is $0 \le \hat{n} \le 2 - 2^{-15}$ = 1.99996948242188 .

the relationship between the two is:
$$ \hat{n} = n \cdot 2^{-15} $$
or
$$ n = \hat{n} \cdot 2^{15} $$ .

note that the "15"(or "-15") in the exponents is the same as the "15" in the "Q1.15" notation.  the scaling factor is directly related to the number of implied bits to the right of the binary point.

if, instead it was a 16-bit signed, we might use "s.15" or "Qs.15" notation to refer to the format, the range of the signed integer $n$ is $-2^{15} \le n \le +2^{15} - 1$.  and the range of the implied fixed-point value, $\hat{n}$, is $-1 \le \hat{n} \le 1 - 2^{-15}$ = 0.99996948242188 .

and the mapping between the two is the same as above:
$$ \hat{n} = n \cdot 2^{-15} $$
or
$$ n = \hat{n} \cdot 2^{15} $$ .

the only difference is the use of unsigned int vs. int in your C code.  they both have the same number of bits to the right of the binary point and they both have the same mapping.

now when C does addition or subtraction or multiplication or division of these values, it is operating on the integer representation of those values.  but the integer representation is the same as your implied value, except for the scaling factors.

so if you want to divide some numerator $\hat{n}$ by some denominator $\hat{d}$ to result in some quotient $\hat{q}$, you must first establish what the ranges of all three numbers are.  then establish how many bits are to the right of the binary point for each number.  then map those numbers $\hat{n}$ and $\hat{d}$ to their integer values $n$ and $d$, and determine what the range of the integer quotient is before it is adjusted to represent $\hat{q}$.  there will be a scaling factor that will pop out at you (a power of 2) which will imply shifting of either $n$ or $d$ (most likely $n$) before the integer division is done.


Different representations of images and their purpose I'm doing a small overview of different image representations used in image processing and computer vision. For every representation, I would like to have (at least) the following information:

the name of the family of representations
some basic characteristics of the representation family (what information is retained in the representation)
what tasks are they typically used for
what are different concrete types of representation from this family
references for different representations and references for their usages

In addition to the representations I already know / have information about it (as an answer listed below), I am mostly interested in compressed domain (e.g. Fourier domain, wavelets), but knowing about other representations, e.g. edge-based, or whatever else I don't know about yet would be great.

If at all possible, I would appreciate if you could provide references (e.g. the article that proposed the representation, or a book explaining the representation), as indicated with the reference-request tag.

Disclamer

I know list-questions are usually discouraged, but we had something similar about segmentation and thresholding, and it worked well.

If the info I already provided is more suited as an answer, I could edit it out and post it as an answer. Also, I'm thinking maybe this would be a good question for community wiki.

These are some representations that I already gathered information on, so I'll list them here.

Pixel based representation

simplest to define, simple neighborhood relations between elements (4,8,6-connectivity (1))
only local information for each element
big number of elements in the representation (used for displaying the image, applications with a lot of noise/details: e.g. medical imaging)

Block-based representations

image divided in a set of (rectangular) arrays
number of elements slightly smaller than with pixel-based, still only local information
both gray-scale and binary images
uses in compression ((2), (3), (4)), segmentation ((4), (5)), extracting different image features ((5), (6), (7))

Region based (sometime also super-pixels) representations

regions not rectangular; formed by grouping similar and connected pixels, usually using over-segmentation
adjacency information between regions (usually as RAG (region-adjacency graph (8)) or combinatorial map (9))
normalized cuts (10), graph-based segmentation (e.g. by Felzenszwalb and Huttenlocher (11)), watershed-based ((12), (13))
number of elements in the representation reduced while accuracy is kept
used for object detection and segmentation, but different unions of multiple regions have to be considered

Hierarchical representations

propose most likely unions of regions of region-based representations
image representation at different scales
only regions present in the representation have to be examined, can be constricted by scale
min-/max-tree ((14), (15)), tree of shapes ((16), (17)), binary partition tree (18), $\alpha$-tree (19), ...
used for: object detection (20), video segmentation (21), image segmentation and filtering ((22), (23)), image simplification (24)


Behaviour of System by Impulse Response In a book by Alan Oppenheim, it is given that for a LTI system, its characteristics are completely determined by its impulse response. But if impulse response of a system is basically an input output relation and impulse response can be different for the same system depending upon our choice of output, so how does it completely characterize the system? Isn't it  contradictory? And is complete characterization of system also include stability of system?

The impulse response characterizes an LTI system completely, in the sense that the response to any input can be computed from the input and from the impulse response  by convolution. That's a consequence of the system being linear and time-invariant.
If we define stability as bounded-input bounded-output (BIBO) stability then stability can be determined from the impulse response. For the system to be stable, the impulse response must be absolutely integrable.
Note that when describing systems by impulse responses, or, equivalently, frequency responses or transfer functions, we focus on the input-output relation. We do not consider the internal dynamical behavior of the system. A state-space model is needed if we not only want to describe the system's input-output behavior but also its internal behavior.


THD requirement for generated sine I need to generate a sine signal using a DAC with a certain THD specification, say 0.1% @ 1KHz. As a first approach I’m considering it to be an ideal DAC along with an ideal LPF for aliasing removal in order to define the minimal resolution in bits so I could achieve the desired THD. To optimize the quality consider I’ll generate the signal using the full scale. Is there a direct method to find the optimal resolution?

I really appreciate any insights!

This can't be the way to go: If you have an ideal LPF, one-bit DAC is always sufficient. (Hint: Draw the spectrum of a square wave. Draw the pass band of your LPF on top.)

You don't have an ideal LPF; or you might want to produce a range of sine frequencies (not only a single one) spanning more than an octave, so that you can't just use a fixed LPF + 1 bit DAC.

The easy way to calculate THD is realizing that it's a relation between powers. So, you need to know how much "error power" there is in a quantized signal. There's this term quantization noise that captures that – and it's well-known for sinusoidal signals!


Ripple in a SciPy Butterworth SOS Filter I have a script using SciPy for checking the RMS of various Butterworth Bandpass Filters of varying orders.
I would expect the RMS values to increase and decrease linearly and consistently as you get closer and further from the passband, but this is not what happens at all. The RMS fluctuates dependent on the frequency, for instance, in the example below, 162Hz has a much higher RMS than 158, despite 162 being further from the high-pass cut-off of 133Hz. This appears to be a cyclical pattern, and independent of the order, but I am not good enough with matplotlib to create a pretty chart that visually represents this issue. The exact cycle also appears to be very dependent on the buffer size.
I am using https://www.szynalski.com/tone-generator/ to test tones, and here is my code:
import numpy as np
import scipy.signal
import pyaudio

def normalize(block):
    array = np.frombuffer(block, dtype=np.float32)
    return array

def get_rms(samples: np.array) -> float:
    """Get the RMS of an array of audio samples

    Args:
        samples: the samples to get the RMS from

    Returns:
        float: the RMS
    """
    samples_array = np.array(samples)
    return np.sqrt(np.mean(samples_array ** 2))  # type:ignore

def design_filter(lowcut, highcut, fs, order=3):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    sos = scipy.signal.butter(order, [low, high], btype="band", output="sos")
    return sos

def main():
    sample_rate = 44100
    buffer_size = 2048

    filters = {}
    for i in range(10):
        sos = design_filter(101, 133, sample_rate, i)
        zi = scipy.signal.sosfilt_zi(sos)
        filters[i] = [sos, zi]

    stream = pyaudio.PyAudio().open(
        format=pyaudio.paFloat32,
        channels=1,
        rate=sample_rate,
        input=True,
        frames_per_buffer=buffer_size,
    )

    update_every = 3
    update = 0
    while True:
        block = stream.read(buffer_size)

        if update_every == update:
            update = 0
            samples = normalize(block)

            # blank out terminal
            print(chr(27) + "[2J")
            # move cursor to to left
            print(chr(27) + "[1;1f")

            print("rms")
            print(int(get_rms(samples) * 200) * "-")

            for order, tup in filters.items():
                print(f"Order: {order}")
                bandpass_samples, zi = scipy.signal.sosfilt(tup[0], samples, zi=tup[1])
                tup[1] = zi.copy()
                print(int(get_rms(bandpass_samples) * 200) * "-")
        else:
            update += 1

if __name__ == "__main__":
    main()

UPDATE
After looking at some more examples and messing with my code, it appears that it is more related to the buffer size than I initially thought.
If buffer_size == sample_rate the ripple (?, oscillation? noise?) does not happen at all (this is wrong, see below). Also, there is seemingly a relationship between factors of buffer_size and sample_rate that relates to the division. Not sure why this is, and I am able to work around it, but if anyone knows why this happens it would be good to know, for sure.
UPDATE 2
It appears that the above is incorrect. There is some transient noise that appears every period. Having a larger period decreases the relative impact of that noise, but it does not eliminate it, and it does not allow for a work around of the issue by increasing the buffer size and then chopping up the filter output later to get the granularity back.

I can find 2 errors at a quick glance :
1 - You should discard the transient when measuring the steady-state gain at for a given frequency $f$. I haven't checked how long the transient last, but you should factor it in your RMS measurement.
2 - The RMS calculation only works when the number of samples corresponds to a whole number of periods. Otherwise, the calculation is biased.
Edit :
Period in this context means the signal period. As you are in the discrete domain, this is the ratio between the sampling frequency and the signal frequency
$T = \frac{f_s}{f}$
For transient, see this site
https://www.dsprelated.com/freebooks/filters/Transient_Response_Steady_State.html
Edit 2 :
Assuming your signal frequency is 1 kHz
$T = \frac{44.1 kHz}{1 kHz} = 44.1$
Since you have 44.1 samples per period, the RMS calculation with only 44 samples will have some bias. However, if you use 441 samples, you will have exactly 10 signal periods and then the RMS calculation will have no bias.
As for the transient, you should plot the time response. You should be able to evaluate how long the transient lasts.


Why Filter Kernel is flipped during convolution process? Why Filter Kernel is flipped during convolution process? What are the consequences if the kernel is not flipped?

Note: Filter Kernel means: Impulse response or convolution kernel.

In addition to the linked question, here is a good one from MathWorks:

The time to flip in so miniscule it's inconsequential. If you go through the theory (linear systems theory) you'll understand.
If you want to understand the causality, imagine a filter that's a right triangle, with the small pointy end on the left and the steep edge on the right. Now shift that past a delta function. If you don't flip and just start shifting your filter past it will hit the tall steep edge first and then ramp down. So then your response is opposite to your filter, but they must be the same if it's a delta function.
Basically it's because time goes along the x axis with the small time values on the left and the big (later) time values on the right. So if you start shifting in, you're having the big time values hit your signal first, which is not right (causal). So you have to flip it to make the small time values shift in first.


Are There Good Methods for Converting Color Image to Gray Scale Image? Converting color image to gray-scale image is often the first step of many image processing tasks. There exists many gray-scale image conversion methods. Most them, however, is to give R G B bands a weighting factor and then the gray-scale image is the weighted sum of R G B bands. I think this is not right in many cases. For example, the most straightforward of way of conversion is to average the R G B bands: gray=0.33*(R+G+B). If we use this formula for generating gray-scale image, then a red object with pixel value R=250 G=0 B=0 will have the same value with a green object (R=0,G=255,B=0). Any ideas of smart color conversion? Thanks.

There are method which are called Contrast Preserving Decolorization.
Those methods are built to keep the contrast in the Color image in the converted image.
They are based on Optimization Problem Solving which creates a Contrast Model in the color space and in the gray scale space.

For instance, have a look on that:


What is the relation between noise (AWGN) variance and number of antennas? I'm working in two Matlab scripts trying to simulate a MIMO system using OFDM (downlink) and SC-FDE (uplink).

My question is: for different numbers of transmit and receive antennas, what is the relationship between the AWGN variance and the number of antennas?
Does it depend if it is uplink or downlink?

I'm using matricial equations like: Y = HX + N , with matrix dimensions: Y(NR x length) , H(NR x NT) , X(NT x length) , N(NR x length)

Usually, the model is that every receiver antenna picks up noise of the same PSD, $N_0/2$. This means that the variance of noise samples in each quadrature branch is $N_0/2$. This noise is local to each antenna and is independent of the rest of the system.

To calculate the total amount of noise in the receiver, you multiply by the number of receiver antennas. This is usually not needed, though, since the figure of merit is the ratio of the energy per information bit (or symbol) and the variance of the noise affecting that bit (or symbol).


Beginner Attempting FFT Signal Filtering With Numpy I've tried looking around for information on this, but I'm really out of my league here. I'm a guy who likes to fool around with Python, and I wanted to make a program that would filter an audio file. I'm using Python and NumPy, with the scipy.io.wavfile module for importing and exporting Wave files. I've gotten it messing around with volume, but not filtering. Here's what I have so far.

rate, data = wavfile.read('./TriLeftChannel.wav')

filtereddata = numpy.fft.rfft(data) # FFT Filtered data
freqdata = numpy.fft.fftfreq(data.size) #Frequency data
filtereddata = AudioFunctions.Filter(filtereddata, freqdata, data, rate) # Filter the data

def Filter(filtereddata, freqdata, data, rate):

    #fftchunks = (rate / 2)# + 1

    x = freqdata[:len(data) / 2]

    for f in range(len(x)):

        if x[f] > 0.1:
            filtereddata[f] = 0.0
            filtereddata[len(data) / 2 + f] = 0.0

    return filtereddata

In that function, filtereddata is the FFT'd data, freqdata is the frequency data that I got with fftfreq(), and data is the wave file itself, 'bare'. Rate is the sampling rate (though I don't use it). This function doesn't actually filter the frequencies (although I know it's a hard filter and no filter should really be this harsh). After I'm done, I output the file with

filteredwrite = numpy.fft.irfft(filtereddata)

filteredwrite = numpy.round(filteredwrite).astype('int16') # Round off the numbers, and get ready to save it as 16-bit depth file (?)

wavfile.write('TestFiltered.wav', rate, filteredwrite)

Since I'm kind of struggling even just to have gotten this far, I was wondering if anyone could give me any pointers or a kind of beginners' tutorial for FFT? Of course, any help would be greatly appreciated.

EDIT: Fixed indentation issues.

Two potential problems with your approach:

You are computing the FFT on your whole signal, which will be terribly inefficient if your input data gets too large, and prevents you from using windowing. If you want to do frequency-domain modifications of your signal, consider using a Short-Term Fourier Transform, processing the resulting FFT frames, and resynthesizing a time-domain signal using overlap-add.
FFT -> zeroing FFT coefficients -> IFFT, especially without windowing of the input data, is seldom used for filtering as it will yield a filter with many unwanted characteristics (side-lobes + non-causal).

A potential problem in your code:

I don't understand why you zero the filtereddata[len(data) / 2 + f] coefficients. I assume it's intended to zero the "negative frequencies" but then the array index is wrong, and you are using numpy.fft.rfft which has a different output layout from numpy.fft.fft (FFT coefficients for negative frequencies are not returned). I'm surprised this runs at all, your code should actually cause an index out of bounds error.

My suggestion would be for you to learn about FIR / IIR filters, design a suitable filter matching your constraints, and apply it to your data. There is a well-established body of work, software tools and practices on FIR / IIR filters and their design. Once you have the coefficients, for small sequences and/or if complexity is not an issue, there's scipy.signal.lfilter to do the job. As suggested by Jason R, applying a long FIR can be made very efficient computationally with FFTs + overlap-add.


FIR filter (from zeros to filter coefficients) Assume $F=1$ $(T=1)$. A 2nd order FIR filter with real coefficients has a zero in $H(e^{j\omega T})$ at $z_0=e^{j\omega /6}$. In addition, $|H(e^{j\pi T})| = 1$.

What are the coefficients of the FIR filter?

current progress

I'm having trouble understanding the $H(e^{j\omega T})$ function.
Tried to start from

$$y[n] = \sum_{i=0}^2 b_i z[n-2].$$

Does that fact that it has a zero in a point means I must apply zero-pole rules here?

First of all a discrete filter having a zero somewhere means it's Z-transform (which is a polynomial) acquires the value 0 there.

If you remember, there is a theorem in calculus saying: for a polynomial to have real coefficients any complex roots which are not on the real line must come in complex conjugate pairs. So if $z_0 = e^{j\omega /6}$ then $z_1 = e^{-j\omega / 6}$ and that gives : $$(z - z_0)(z - z_1)C = (z - e^{j\omega /6})(z - e^{-j\omega /6})C$$

Then to find the constant $C$ solve the equation with the absolute value of the polynomial equals 1 for $z = e^{j\pi}$.


Extract sharpness from FFT I've got in total 6 images, from 2 cameras (thus 3 per camera) storing 16 bit data.

What I need to feedback to the user is like some sort of sharpness coefficient, preferably live while he is adjusting the focus of the whole system (since its looking into the eyes, it depends on the dioptres of the patient).

What I am currently using is a contrast coefficient combined with standard deviation, which should work fine, but with a change of the focus the field of view changes slightly, which can impact the contrast.

Now I wonder when I move to use the FFTW library to get the FFT of my image, it results in an array of fftw_complex (which is double[2]).

For my purpose, it should be fine to calculate the average of the resulting array, but I don't understand how to mangle real / complex components in this. Do I iterate over the whole array, calculate sum equalling sum of double[0] and double[1], or am I totally wrong here? I don't really understand why there is a real and a complex component. I understand that the FFT results in coefficients for different frequency sinus waves, but the frequency of the wave depends on its position in the array, while the amplitude should be the real component. Is the combination of real and complex variable responsible for the direction in the image, or what is it used for, and how would I use it to get an idea of an "average" frequency (sharpness) in the image?

As you said, in the fourier transform, the frequency of the wave depends on its position in the array. But, the amplitude of the wave is not always the real component, it is actually the magnitude of the complex number ($\sqrt{real^2+imaginary^2}$).

If your initial image is perfectly symmetric along the exact center, you can get real coefficients, but it's something which seldom appears in image analysis.

Edit: Here's more details... Couldn't find the picture I was refering to...

The discrete Fourier transform, the basis of every FFT algorithm, is defined as such:
$$X_k = \sum_{n=0}^{N-1}x_n\cdot e^{-i2\pi kn/N}$$

From Euler's formula, a complex exponential can be separated in a complex sum of cosinus and sinus:
$$e^{ix}=\cos(x)+i\sin(x)$$

If we now look back at the Fourier transform, it can be separated into two sums:
$$X_k = \sum_{n=0}^{N-1}x_n\cos(-2\pi kn/N) + i \sum_{n=0}^{N-1}x_n\sin(-2\pi kn/N)$$
As a sidenote here, if your image (or signal) is perfectly symemtric about it's center, the cosinus will resonate strongly with your image while the sinus will not output anything. It will result a completly real image. It's the other way around for a perfectly anti-symmetric signal where one half is the negative of the other.

Now knowing this, it's straightforward to compute the amplitude of the resulting FFT, using the formula I gave in the first paragraph.

To get the frequency of the image itself is less obvious. We would expect the result to look somewhat like this, with the zero-frequency component at the center. It is, however, not the case.

The "standard" FFT algorithm will put the zero-frequency component in a corner (I can never remember which one, top-left or bottom-left!), such that the array indice [0,0] corresponds to $k_x=k_y=0$. But wait, once you cross N/2, the pixels represents negative frequencies, counting backwards. A good reminder is Python's scipy.fftpack.fftfreq() function found fftfreq().

My personal tips to take a quick look at an image's spectrum is to shift it then display the logarithm of it. When it's not shifted, it's hard to see asymetries. When looking at the non-log version, it's usually just a big black square with a small fuzzy grey dot...


Periodicity and fundamental frequency detection algorithms I have some time series datasets which are periodic, and I'd like to pick out the fundamental frequency of the data and take the average over a fundamental period.

Questions:

1) Would something simple like taking the FFT of the time series data and looking for the strongest peak frequency do the trick, or are there more complicated period detection algorithms out there? For what it's worth, I'm trying to do all this in MATLAB.

2) There are some edge effects, particularly at the beginning of each data set. Would I have to take this into account?

This is essentially the problem of Pitch Detection.  A fairly simple approach is to compute the autocorrelation (FFTs can be used to do this quickly),  and look for the peak with T > T_min


Channel impulse reponse multipaths I find this image :

Does it mean that at 0ns there is one path, 10ns second path, ... and approximately at 130 ns, there is a power path?

In the multi-path scenario we have occurrence energy from (typically) all directions, and a channel could be a Rayleigh fading channel; no LOS component or Rician fading channel; a strong LOS component is present.

The figure you referred to contains a dominant component with a delay spread of 130 ns and with responses from other paths occurring with different delay spreads. Also, the key to the figure mentions an RMS delay spread of 65.9 ns, which can be justified by noticing the higher power distribution between 0 to 100 ns delay spread.

Take a look at the picture below to better understand the CIR for a multipath scenario. In this case, the CIRs from multiple antennas will be combined to form the CIR of the antenna array.


Phase distortion of frequency warping by using Praat I am trying to do speaker conversion by frequency warping method. I modify the function Spectrum_shiftFrequencies() in Praat by scaling the FFT spectrum of a voice. When one scaling factor is multiplied to spectrum, the timbre of source voice is changed and the voice quality is good. However, using multiple scaling induces several frequency bands of echos into the synthesis audio. I guess that the multiple scaling method breaks the phase information of a spectrum.
Fig.1 single scaling and multiple scaling

Here is the code that I modify in Spectrum_shiftFrequencies() in Praat.

for (long i = 1; i <= thy nx; i++) {  //scan target frequency component
    double thyf = thy x1 + (i - 1) * thy dx;  // target freq
    double myf =  0;                 //  from source freq

    if(thyf>0 && thyf <=target_formant[0]) // 2 scaling
        myf=thyf*source_formant[0]/target_formant[0];
    else
        myf=source_formant[0] + (thyf-target_formant[0])*(my xmax-source_formant[0])/(my xmax-target_formant[0]);

        if (myf >= my xmin && myf <= my xmax) {
            double index = Sampled_xToIndex (me, myf); // index:  source freq index

            thy z[1][i] = NUM_interpolate_sinc (my z[1], my nx, index, interpolationDepth);
            thy z[2][i] = NUM_interpolate_sinc (my z[2], my nx, index, interpolationDepth);
      }

}

The real and image parts of target spectrum are interpolated from source spectrum. After checking the amplitude spectrum of synthesis audio, the contour of spectrum is changed as the expected scaling function. The phase spectrum is also changed, but that induces echo.

Does the phase spectrum require special process in this multiple scaling (nonlinear) conversion? or How to eliminate the echo ? or What is the correct phase conversion of this frequency mapping ?

Any help is appreciated. Thank you.

After research works, I find that the answer is the pitch shift. The terminology is quite confusing. The proper term should be "frequency shift" because it changes the spectrum of signal, not only pitch but also formant.

Some useful resource is listed bellow:

Pitch Shifting Using The Fourier Transform

guitar pitch shifter


Multi-scale in Object detection? I want to ask about MultiScale Object Detection. For the multiscale, if we have an original image do I need to zoom-in at different scales in the image at a particular point or zoom-out at different scales in the image at a particular point?

Multiscale object detection generally consists in smoothing the image with larger and larger kernels, to check  which objects persist across different scales. So, zoom-out, mostly. You now can have access to the scale-space representation of the image, as shown below.

One  reason is that, without more knowledge, one does not zoom in beyond the pixel width. As illustrated (with humour) in the following picture.

Meanwhile, in certain conditions, subpixel resolution or superresolution is possible.


Ultrasonic data transfer I am trying to send a short amount of data (10-15 characters) from one mobile device to another via (inaudible) sound within 1-2 seconds. I'd like to be able to deal with normal external noises (chatter, traffic, etc).

I have a software engineering background, but no experience working with sound or audio. Is there a known/documented solution to address this problem of reasonably reliable ultrasonic data transfer? If you know of a guide or some document that'll explain the concepts to me as I develop out a solution that would be really awesome.

I've researched a little bit on this topic, but have been unable to find something concrete for ultrasonic communication. Any direction would be much appreciated (and yes, I know I'm a noob here). Thanks.

Even using the highest audio frequencies that a typical mobile device can generate (which are below the ultrasonic range), there is no guarantee that "normal external noise" won't interfere.  Transient impact noises can be quite wide-band.  Therefore one may have to try various forms and levels of redundancy or Forward Error Correction coding (Hamming, Reed-Solomon, etc.), until the reliability probabilities reach a suitable level for your application.  For many (most?) modulation schemes, just lowering the bandwidth sent/used can also improve the reliability probabilities of the channel.


Discrete Time Convolution Evaluation I'm trying to solve a problem on convolution from Alan V.Oppenheim:
Find the convolution output $y[n]$ for the following signals:
$$x[n]= u[n]\quad\text{and}\quad h[n]=a^{n}u[-n-1], \ a>1 $$
I started the evaluation:
$$y[n]=\sum_{k=-\infty}^{+\infty} u[k]a^{n-k}u[-n+k-1]$$
considering that $u[k]=1$ for $k\ge0$ and $u[-n+k-1]=1$ for $k\ge n+1$
which I evaluated to
$$y[n]=a^{n}\sum_{k=m}^{+\infty} a^k$$ where $m=n+1$ could be $<0$ or $>0$
and I tried to evaluate for $m>0$ which is the same as $n>-1$, which evaluated as:
\begin{align}
y[n]&=a^{n}\sum_{k=m}^{+\infty} a^m\\
&=a^{n}\left[\left(\sum_{k=0}^{+\infty} a^m\right)-\left(\sum_{k=0}^{m-1} a^m\right)\right]\\
&=a^{n}\left[\left(\frac{1}{1-a^{-1}}\right)-\left(\frac{1-a^{-m}}{1-a^{-1}}\right)\right]\\
&=\frac{a^{n-m}}{1-a^{-1}}\\
&=\frac{a^{-1}}{1-a^{-1}}\tag{since $n-m=-1$}
\end{align}
but when I evaluated for $m<0$  which is $n\le -1$ I am facing a problem:
$$y[n]= a^{n}\left[\left(\sum_{k=m}^{-1} a^m\right)+\left(\sum_{k=0}^{+\infty} a^m\right)\right]$$
How do I evaluate the first summation? I mean am I to consider $k=-m$ since $m<0$?

The multiplication of the two unit step sequences $u[k]\cdot u[-n+k-1]$ is only non-zero if both sequences are non-zero. This means that the condition $k\ge 0$ as well as the condition $k\ge n+1$ must be satisfied. So you have two cases: for $n<=-1$ you have to evaluate the sum with the lower limit $k=0$, and for $n>-1$, you have to evaluate the sum with a lower limit $k=n+1$.


Discrete approximation to Gaussian filter One method of producing a discrete approximation kernel to a Gaussian filter of variance $\sigma^2$ is to assume a cutoff of $5\sigma$. This seems to suggest the following procedure:

Create a uniform grid $\mathbb{x}$ of size $11(=2*5+1)$ between $[-5\sigma,5\sigma]$
Evaluate $\mathcal{N}(\mathbb{x};0,\sigma)$ and obtain the discrete approximation.

The blue lines in the second plot below depict the discretized values(for $\sigma=1.5$).
In this case, the result is not a proper filter kernel since the sum of values is not $1$.

But MATLAB's fspecial does things a bit differently. For the call fspecial('gaussian',[1 11],sigma)

Create a uniform grid $\mathbb{x}$ of size $11(=2*5+1)$ between $[-5,5]$
Compute $\mathbb{x}'= e^{\frac{-\mathbb{x}^2}{2\sigma^2}}$
Normalize $\mathbb{x}'$ so that it sums to $1$

The green lines in the first plot below the values as computed by MATLAB's fspecial(for $\sigma=1.5$).

What is the difference between both methods ?

The difference between the two should be slight but the MATLAB one has the property that it preserves the signal energy content. That is a property that you expect to have when using a Gaussian filter.


Understanding a Bainter Circuit (chebychev filter) As part of a class in Signal Processing, I am building a 3rd order Chebychev band reject filter. We implement this by using three cascading Bainter circuits. Although not part of the class, I have a question about the gain of the Bainter circuit.

I am trying to write a script which will automate the component selection using the corner frequencies and maximum overall gain as design rules, but have some issue with calculating the overall gain.

To calculate the overall gain of a Bainter stage, would I simply work out the individual gains of the three op-amp sections? The overall gain would then be the product of the three individual gains?

To calculating the overall gain of a Bainter stage, would I simply
  work out the individual gains of the three op-amp sections. The
  overall gain would then be the product of the three individual gains?

The short answer is: Yes, you can (probably) analyze them individually.

When asking what happens when you cascade multiple analog filter stages, the questions to ask are: what is the source impedance of the first stage, and what is the load impedance of the second stage?  If a circuit stage has a large and complicated output impedance, then loading it with another stage can modify its behavior.  When working with passive filters, this is a big problem: unless the load impedance of each stage is significantly greater than the source impedance of the prior stage, cascading passive filter sections will result in complicated changes to the behavior of each stage.

One of the attractions of op-amp based circuits is that op-amps generally have a very low output impedance; for the ideal op-amp, there is zero output impedance.  Moreover, the op-amp inputs themselves generally have very high input impedance, ideally infinite.  This means that circuit sections whose outputs are op-amp driven can generally be cascaded without having one stage change the behavior of another.

Consider this schematic of a Bainter notch (taken from an Analog Devices publication):

The "notch out" is driven by the output of an op-amp.  Thus this circuit will have very small output impedance.  In other words, the voltage at "notch out" will be relatively insensitive to the load that is connected.   This output impedance will almost certainly be much lower than the input impedance.

Thus, in the design phase, you can analyze several cascaded notch circuits separately, and simply multiple their transfer functions together.  After producing a design in such a way, you might want to simulate the entire circuit in SPICE to check for behaviors due to op-amp nonidealities, etc.

References

"Active filter has stable notch, and response can be regulated", in Electronics, October 2, 1975.  Original article describing the circuit.
"How to Write Node Equations for Active Circuits", by Mr. Bainter himself.


Blurring Non-Uniform Areas I have an image with arbitrary area (the boundary is shown in gray):

I would like to blur just the image within boundaries.

Ordinary blur causes the background proliferating to image and vice versa. This results in dark areas near boundaries of the image:

One simple solution would be to create a binary (0/255) mask of the area you want to blur. Then blur the source image BITWISE AND mask, blur the mask with the same filter and divide them. As pseudocode:

(filter(source & mask) / filter(mask)) & mask


What criteria should be used for smoothing Actually I have some measurement which I want to get rid of noise
I want to use different filter techniques but I am wondering what criteria I should use to check if I am removing noise or not
or if my selected parameter is optimal.

My measurement is mass spectrometry data.

You want to look at the innovation sequence; the difference between the filtered value (the estimate or forecast) and the measured value. If you have a perfect filter that removes the noise, the innovation sequence should be white noise (Dirac impulse autocorrelation).

https://en.wikipedia.org/wiki/Innovation_(signal_processing)

Here is some MATLAB code that illustrate these properties for some white sequences.

N = 1e5;
n1 = -1 + 2*rand(1,N); % uniform noise
n2 = n1 + -1 + 2*rand(1,N); % triangular noise
n3 = randn(1,N); % Gaussian noise

[c1,l1] = xcorr(n1); [c2,l2] = xcorr(n2); [c3,l3] = xcorr(n3);
figure(1), plot(l1,c1,l2,c2,l3,c3) # autocorrelation (Dirac delta)
figure(2), clf, hold on
pwelch(n1,kaiser(N/1000,10),N/2000) # power spectral density (constant)
pwelch(n2,kaiser(N/1000,10),N/2000)
pwelch(n3,kaiser(N/1000,10),N/2000)


NDA timing estimation for mutli-level QAM(RC pulse) What non-data0aided techniques are used to estimate the timing for multi-level QAM? For single level QAM(so 4 QAM) you can use variance minimization or power maximization but I'm not sure these techniques carry over to multi-level transmission models. For example the variance minimization technique would need substantial modifications to handle the different power levels(even more so when you consider lack of phase synchronization).

How is this addressed?

Consider using a Gardner Timing Error Detector which in the following form is usable for higher order QAM:

$$TED = I_n(I_{n+1}-I_{n-1}) + Q_n(Q_{n+1}-Q_{n-1})$$

Where $I_{n-1}$,$Q_{n-1}$,  $I_{n}$, $Q_{n}$ and $I_{n+1}$, $Q_{n+1}$ are the early, prompt and late QAM samples at 2 samples per symbol. The TED will drive the prompt sample to zero error at the zero crossings of the waveform.

This is from deriving a timing error by taking the real of the complex conjugate multiplication of the signal y[n] with its derivative, a form which comes from simplifying the maximum likelihood timing estimator. The difference $y_{n+1}-y_{n-1}$ is an approximation of the derivative of y[n]:

$$real\bigl(y^*_n (y_{n+1}-y_{n-1})\bigr)$$

This error term would typically be the input into a loop filter in a timing recovery loop, and the output of the loop filter could select the appropriately timed output of a polyphase resampling structure such that the sampling rate is never actually changed (2 samples per symbol throughout).

Assuming the RC pulse is root-raised cosine (RRC) in the transmitter, and then again in the receiver, the Gardner TED has higher performance using the samples prior to the RRC filter in the receiver: this makes sense as the TED derives the timing error or S-curve based on the symbol transitions- if you observe the zero-crossing jitter in an eye diagram of RRC pulses compared to RC pulses with the same pulse-shaping parameter, you will see that the RRC has less jitter. The second RRC filter adds more jitter at the benefit of eliminating the inter-symbol interference at the symbol sampling locations.

That said, the waveform can be prefiltered prior to the timing error detection to minimize zero crossing jitter (this filtered waveform is used only for timing and is not the received waveform for data recovery due to ehanced ISI at the correct symbol decision locations) Here is an interesting paper on the use of such noise minimization pre-filters for QAM with the Gardner Ted:

Aldo Nunzio D’Andrea, Marco Luise, Optimization of Symbol Timing Recovery for QAM Data Demodulators, IEEE Trans on Communications, Vol 44 No 3 March 1996.

https://ieeexplore.ieee.org/document/486334

Further, here is an interesting graph I created showing how robust the TED is without carrier/phase synchronization. This shows how the TED performs as an SNR metric with carrier offsets normalized to the symbol rate, showing only 3 dB degradation with carrier offsets as much as a quarter of the symbol rate. (The "Gain" is from the slope of the discriminator, while the noise is the pattern noise of the discriminator. For the following two graphs these were just consistent metrics to compare various timing discriminators such as here the Gardner TED versus the M&M synchronizer.)

For comparison here is the same plot with the Mueller and Mueller (M&M) synchronizer. Notice the $10^{-4}$ scaling on the horizontal axis! This is overlaid on the plot above with the red line showing the same TED curve, just zoomed in on the horizontal axis; What this shows is that the M&M is much more sensitive to carrier offset (as expected since it operates on decisions) but has lower noise once the carrier is synchronized.

Here is a simplified plot showing how the Gardner TED estalishes a timing discriminator "S-curve" from the product of the signal with its derivative. This shows how the Gardner uses the zero crossing transitions to establish the discriminator curve (the actual curve is sinusoidal for raised cosine waveforms) and importantly how the lock point is at the zero crossings (blue samples) while the decision samples would be the other samples shown in red (with 2 samples per symbol). For higher order QAM the functional operation is the same with substantially higher pattern noise. (The timing loop BW would need to be reviewed to confirm the resulting discriminator noise is sufficiently averaged to not be an issue).


How does Songify work? Does anyone know how does the "Songify" app algorithm work?

There are 2 inputs:

speech (or amateur singing), and
underlying music.

Somehow the speech is altered and mixed into the music. Which DSP algorithms are needed to achieve this effect?

Pitch-shifter?
Time-stretcher?
Fundamental (F0) frequency detection?
Auto-tune?
Else?

What would be the steps?

Songify seems use Prosodic Modification (Pitch, Time) for monophonic signals:

First you need get the pitch contour and pitch duration from musics you want to follow
(use one pitch track for monophonic signal like YIN,
AMDF, Auto-Correlation)
Split your imput signal in voiced/unvoiced
Apply Pitch Scale Modification in voiced parts to match your extracted pitch contour (you can use PSOLA technique here)
Apply Time stretch to match your extracted pitch duration (you can use PSOLA technique here)


Hardware requirements for ambient active noise cancallation system using adaptive LMS Background

I'm working on a hobby project to reduce ambient low-frequency noise generated from HVAC equipment within a small area of a room. The noise is predictable and concentrated between 55Hz and 120Hz. Aggressive passive noise attenuation efforts have failed at these frequencies, and so ANC may be the only option. I have basic analog mixing equipment, microphones and speakers, and am now researching DSP and related hardware requirements for ANC.

Algorithm

My preliminary research has led me to believe that a Filtered Least Mean Squared (FxLMS) algorithm may be the best choice for an application like this. Because it would be implemented in a "open" space, my understanding is that the secondary path filter of the FxLMS (as opposed to a simple LMS) would be critical. The related diagram would be:

Where:

x(n) is the reference signal recorded via a mic near the HVAC noise source, fed through a preamp and available as a line-level signal.
e(n) is the error signal recorded via a mic in the cancellation area, fed through a preamp and available as a line-level signal.
y(n) is the cancellation signal that will be sent as a line-level signal to a powered sub with an integrated LPF near the cancellation area (filter output).
d(n) is the ambient noise to be cancelled after traveling from the source the the cancellation area.

My understanding is that this is conceptually a pretty straightforward ANC setup, and it is referenced in some detail in numerous publicly available resources (e.g. here and here). However, what is not clear to me is the hardware requirements to implement the DSP filter part (the blue box).

Possibilities

I'm assuming that the real-time synchronization between the sampled x(n)/e(n) signals and the final y(n) output requires a specialized DSP board. I suppose that the algorithm itself could be programmed in Matlab or C (samples are available), but refining, compiling and embedding the result onto a board from scratch seems like something that may require expertise that I do not currently have. So I'm wondering if there are more "black box" solutions available somewhere for this sort of thing.

The closest things I could find are the SigmaDSP Audio Processors, possibly the evaluation versions of this board or this one. If I'm understanding it correctly these come with a huge library of pre-compiled filters that can be chained together and mapped to various I/O via a GUI. One of the available filters seems to be an FxLMS (source), but it's hard to tell which SigmaDSP boards are capable of running that filter or if there are other special considerations that I need to be aware of given an application like this.

So on the hardware side what would be the "lowest barrier to entry" for an application like this? Though I have a background in systems engineering, I do not have extensive experience in audio signal processing. I'm also wondering if this is the kind of problem space that may require significant additional background theory knowledge before tackling.

I would recommend checking out "Design of active noise control system with the TMS320 family" by Sen Kuo [[report link]].

TMS320 chips are specifically designed to perform various ANC tasks. In the report, they even provide code for programming LMS and FxLMS.

In a lot of the literature on real-time ANC, the control system is implemented using these chip sets. So in my opinion, these chips would be a good place to start. Alternatively, I'm sure this system could be implemented with a variety of microcontrollers, but I am not aware of any other hardware specifically designed for ANC.

Side note: I would consider looking into the leaky FxLMS algorithm as well. It's no more complex than regular FxLMS, but provides more stability to the adaptive algorithm (although reduces the steady state performance). The leakage term is typically necessary in real-time practical control systems.


Is there a way to measure each pixel's content in an image somehow? Is there a way to measure each pixel's content in an image somehow?

Like can I interpret each color as some kind of a single value, rather than a triplet?

Perhaps I can think of them as $\mathbb{R}^3$ vectors and then e.g. take their norms?

I would want to start from e.g. making a measure of "color" that's a single real value, rather than a triplet. Then I could perhaps want to make a function that calculates the "change in" color between pixels adjacent to each other. So a derivative of some kind

Why would you need the color to be a single value for that? I'll do an analogy: To measure the distance between to cities, you don't need to map these two cities onto the real numbers; they still have coordinate vectors. All you do is map the difference vector onto a number. And for numbers, it's just the same: It's mathematically impossible to sensibly map them onto a single dimension and still preserve distance. So, you can't.

What you need is just some kind of a norm to map a color difference vector to a "perceptive difference". For many color systems, something like the squared sum of vector components would be a good start.


Proof of DTFT equal to DFT when signal is periodic? I was using the Wikipedia page on the discrete time Fourier transform to understand the connection between DFT and DTFT. The following is claimed in the article - I was wondering if anyone had a proof or a source for the proof for the highlighted part?

Many thanks!

The proof is rather simple if you consider the counter example and I will leave it to the OP to work out the details with these hints that are hopefully helpful.
Understanding the premise of the Fourier Series Expansion is a good starting point: Any single valued analytic function from $0$ to $T$ in time can be decomposed into an infinite series of sinusoids each with a frequency of $n/T$ where $n$ is any integer. Here is the key point:
The result of summing all those components will result in a waveform that matches the original function as expected from $0$ to $T$, but also since each of those components repeat themselves over that interval, then the original function will also repeat over every additional interval of $T$ in time! Consider that any other frequency other than these provided will not repeat as such, then they would disrupt the resulting repetition so therefore cannot exist.
When the beginning of a waveform is not continuous with the end of the waveform, a discontinuity results in the repeated pattern, leading to ringing at the edges in the Fourier reconstruction known as the Gibbs effect. This is consistent with the fundamental requirement that the waveform needs to be single valued and analytic function.


Matlab: Radix 2 FFT with 8 Channels and High Number of Inputs I'm making in Matlab an 8 channel FFT, according to an online tutorial I found at https://www.youtube.com/watch?v=EsJGuI7e_ZQ&t=957s . It works fine for the same number of inputs as shown in the video (N=8), but as I increase N (say for N=128) it no longer works.

I am not sure why it no longer works as I increase N, but I suspect it has something to do with the twiddle factors. I assumed the twiddle factors would have an exponent increasing by 8 after each stage, but maybe they change in a different way as one increases N?

Here is a picture from the YouTube tutorial, at 15:16

Here is a signal flow diagram of my Matlab code

And below is my actual MatLab code. If one uses N=8 (which you can define at top of code) then output is correct. Using N=128 (or other higher values) produces an incorrect result.

I am very confused about what could be causing this problems with the output. Any thoughts?

Much appreciated,

clear all

% Generate input data sequence and plot
N           = 128;
f1          = 10;
num_cycles  = 2;
fs          = f1*N/num_cycles;
x_time      = 0:1/fs:num_cycles/f1-1/fs;
X           = sin(x_time*2*pi*f1);
plot(x_time,X);
title('Input Waveform');

% Split inputs into eight channels
X0_0  =  X(1:8:N);
X0_1  =  X(5:8:N);
X0_2  =  X(3:8:N);
X0_3  =  X(7:8:N);
X0_4  =  X(2:8:N);
X0_5  =  X(6:8:N);
X0_6  =  X(4:8:N);
X0_7  =  X(8:8:N);

% Compute FFT of each channel
X1_0  =  fft(X0_0);
X1_1  =  fft(X0_1);
X1_2  =  fft(X0_2);
X1_3  =  fft(X0_3);
X1_4  =  fft(X0_4);
X1_5  =  fft(X0_5);
X1_6  =  fft(X0_6);
X1_7  =  fft(X0_7);

% Generate Twiddle factors
Wn=exp(-1i*2*pi/N);

% Produce output of first stage of butterfly

for k=0:(N/8)-1
    X2_0(k+1)  =  X1_0(k+1) + (Wn^(k*8)) * X1_1(k+1);
    X2_1(k+1)  =  X1_0(k+1) + (Wn^(k*8+4)) * X1_1(k+1);
    X2_2(k+1)  =  X1_2(k+1) + (Wn^(k*8)) * X1_3(k+1);
    X2_3(k+1)  =  X1_2(k+1) + (Wn^(k*8+4)) * X1_3(k+1);
    X2_4(k+1)  =  X1_4(k+1) + (Wn^(k*8)) * X1_5(k+1);
    X2_5(k+1)  =  X1_4(k+1) + (Wn^(k*8+4)) * X1_5(k+1);
    X2_6(k+1)  =  X1_6(k+1) + (Wn^(k*8)) * X1_7(k+1);
    X2_7(k+1)  =  X1_6(k+1) + (Wn^(k*8+4)) * X1_7(k+1);
end

% Produce output of second stage of butterfly

for k=0:(N/8)-1
    X3_0(k+1)  =  X2_0(k+1) + (Wn^(k*8)) * X2_2(k+1);
    X3_1(k+1)  =  X2_1(k+1) + (Wn^(k*8+2)) * X2_3(k+1);
    X3_2(k+1)  =  X2_0(k+1) + (Wn^(k*8+4)) * X2_2(k+1);
    X3_3(k+1)  =  X2_1(k+1) + (Wn^(k*8+6)) * X2_3(k+1);
    X3_4(k+1)  =  X2_4(k+1) + (Wn^(k*8)) * X2_6(k+1);
    X3_5(k+1)  =  X2_5(k+1) + (Wn^(k*8+2)) * X2_7(k+1);
    X3_6(k+1)  =  X2_4(k+1) + (Wn^(k*8+4)) * X2_6(k+1);
    X3_7(k+1)  =  X2_5(k+1) + (Wn^(k*8+6)) * X2_7(k+1);
end

% Produce output of third stage of butterfly

for k=0:(N/8)-1
    X4_0(k+1)  =  X3_0(k+1) + (Wn^(k*8)) * X3_4(k+1);
    X4_1(k+1)  =  X3_1(k+1) + (Wn^(k*8+1)) * X3_5(k+1);
    X4_2(k+1)  =  X3_2(k+1) + (Wn^(k*8+2)) * X3_6(k+1);
    X4_3(k+1)  =  X3_3(k+1) + (Wn^(k*8+3)) * X3_7(k+1);
    X4_4(k+1)  =  X3_0(k+1) + (Wn^(k*8+4)) * X3_4(k+1);
    X4_5(k+1)  =  X3_1(k+1) + (Wn^(k*8+5)) * X3_5(k+1);
    X4_6(k+1)  =  X3_2(k+1) + (Wn^(k*8+6)) * X3_6(k+1);
    X4_7(k+1)  =  X3_3(k+1) + (Wn^(k*8+7)) * X3_7(k+1);
end

% Merge X4 into final output

c = 1;
for k = 1:N/8
    X5(c:c+7) = [X4_0(k);X4_1(k);X4_2(k);X4_3(k);X4_4(k);X4_5(k);X4_6(k);X4_7(k)];
    c = c+8;
end

% Plot expected FFT and Butterfly FFT
figure
matlab_fft=fft(X);
plot(abs(matlab_fft));
title('Matlab FFT');
figure
plot(abs(X5));
title('Butterfly FFT');

Your problem is probably with twiddle factor adjustment. Below is a working code for 8 channel implementation of FFT (where the last 3 stages are explicitly computed, but not in a butterfly structure). Its derived based on the recursive structure of FFT. Where the splitting algorithm of N-point DFT into 2 N/2 point DFTs is applied for three consequtive stages.

Mathematical definition of the solution is as follows: Assume that $N$-point DFT $X[k]$ is split into two $N/2$ point DFTS according to the even and odd partitioning as usual:
$$ X[k] = X_e[k] + W_N^k X_o[k] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

where the twiddle factor is $W_N = e^{-j\frac{2\pi}{N}}$ and note that the partial DFTs $X_e[k]$ and $X_o[k]$ are of length $N/2$and periodic. Then we can further decompose those $X_e[k]$ and $X_o[k]$ recursively to reach

$$\begin{align}
X_e[k] &= X_{ee}[k] + W_{N/2}^k X_{eo}[k] ~~~,~~~\text{ for } k=0,1,...,N/2-1 \\
X_o[k] &= X_{oe}[k] + W_{N/2}^k X_{oo}[k] ~~~,~~~\text{ for } k=0,1,...,N/2-1 \\
\end{align} $$

And further split those $N/4$ point DFTs into $N/8$ points DFTS which will be the channel signal DFTS.

$$\begin{align}
X_{ee}[k] &= X_{eee}[k] + W_{N/4}^k X_{eeo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
X_{eo}[k] &= X_{eoe}[k] + W_{N/4}^k X_{eoo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
X_{oe}[k] &= X_{oee}[k] + W_{N/4}^k X_{oeo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
X_{oo}[k] &= X_{ooe}[k] + W_{N/4}^k X_{ooo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\
\end{align} $$

Below is an implementation of the above splitting algorithm. Note that your definition of the signals X0_0,X_0_1...,X0_7 is actually correct according to this splitting mode. But the butterfly computations seems not right.

clc;clear all; clear all

% S0 - Generate N point input data sequence x[n] of a sine wave.
% --------------------------------------------------------------
N           = 128;          % Length of signal
M           = 8;            % Number of channels
x           = sin(2*pi*0.123*[0:N-1]);

% S1 - obtain the even-odd partition signals :
% --------------------------------------------
xe = x(1:2:N);              % EVEN part of N-point x[n]
xo = x(2:2:N);              % ODD  part of N-poibt x[n]

xee  = xe(1:2:N/2);         % EVEN part of N/2-point xe[n]
xeo  = xe(2:2:N/2);         % ODD  part of N/2-point xe[n]
xoe  = xo(1:2:N/2);         % EVEN part of N/2-point xo[n]
xoo  = xo(2:2:N/2);         % ODD  part of N/2-point xo[n]

xeee = xee(1:2:N/4);         % EVEN part of N/4-point xee[n]
xeeo = xee(2:2:N/4);         % ODD  part of N/4-point xee[n]
xeoe = xeo(1:2:N/4);         % EVEN part of N/4-point xeo[n]
xeoo = xeo(2:2:N/4);         % ODD  part of N/4-point xeo[n]
xoee = xoe(1:2:N/4);         % EVEN part of N/4-point xoe[n]
xoeo = xoe(2:2:N/4);         % ODD  part of N/4-point xoe[n]
xooe = xoo(1:2:N/4);         % EVEN part of N/4-point xoo[n]
xooo = xoo(2:2:N/4);         % ODD  part of N/4-point xoo[n]

%S1 - Compute 8 of those N/8 point FFTs for the eee-ooo signals :
% ---------------------------------------------------------------
Xeee = fft(xeee,N/8);
Xeeo = fft(xeeo,N/8);
Xeoe = fft(xeoe,N/8);
Xeoo = fft(xeoo,N/8);
Xoee = fft(xoee,N/8);
Xoeo = fft(xoeo,N/8);
Xooe = fft(xooe,N/8);
Xooo = fft(xooo,N/8);

% S1 - obtain N/4 point DFT's from 8 , N/8 point FFTs of eee-ooo signals :
% ------------------------------------------------------------------------
W4 = exp(-1j*2*pi*(4/N));      % twiddle factor for the 3rd stage (Wn/4)
Xee = [Xeee, Xeee]  + (W4.^[0:N/4-1]).*[Xeeo, Xeeo];
Xeo = [Xeoe, Xeoe]  + (W4.^[0:N/4-1]).*[Xeoo, Xeoo];
Xoe = [Xoee, Xoee]  + (W4.^[0:N/4-1]).*[Xoeo, Xoeo];
Xoo = [Xooe, Xooe]  + (W4.^[0:N/4-1]).*[Xooo, Xooo];

% S2 - obtain N/2 point DFT's from N/4 point FFTs of ee-oo signals :
% ------------------------------------------------------------------
W2 = exp(-1j*2*pi*(2/N));      % twiddle factor for the 2rd stage (Wn/2)
Xe = [Xee, Xee] + (W2.^[0:N/2-1]).*[Xeo, Xeo];
Xo = [Xoe, Xoe] + (W2.^[0:N/2-1]).*[Xoo, Xoo];

% S3 - obtain N point final DFT from two N/2 point DFTs of e-o signals :
% ----------------------------------------------------------------------
Wn = exp(-1j*2*pi/N);
X = [Xe, Xe] + (Wn.^[0:N-1]).*[Xo, Xo];

% SX - Display results
% --------------------
figure,stem(abs(X))
figure,stem(abs(fft(x,N)),'g');


how to create synthetic $1/f$ noise? I am writing an app to work with synthetic time series data from a physics experiment. In our experiments we always have $1/f$ noise in our time series, but I haven't been able to find code/packages to generate $1/f$ noise in synthetic time series data. Can anyone tell me how to do this? I have a feeling it should be pretty simple and that I'm missing something obvious, but I don't know what it is. Any advice?

Also, I looked through the list of questions and don't see anything too on point apart from this: making pink noise (1/f) using list of frequencies

this is an old resource, but good.   in fact, i have a very old contribution (another 3 decade old design) that i will repeat here at SE:

A method that Orfanidis (Introduction to Signal Processing mentions came from an old comp.dsp post of mine. (here's a pdf, check Problem B.9)  It's just a simple "pinking" filter to be applied to white noise.  since the rolloff is -3 dB/octave, -6 dB/octave (1st order pole) is too steep and 0 dB/octave is too shallow.

An equiripple approximation to the ideal pinking filter can be realized by
alternating real poles with real zeros.  A simple 3rd order solution that i
obtained is:

$$ H(z) = \frac{(z-q_1)(z-q_2)(z-q_3)}{(z-p_1)(z-p_2)(z-p_3)} $$

where

$p_1$ =  0.99572754

$p_2$ =  0.94790649

$p_3$ =  0.53567505

$q_1$ =  0.98443604

$q_2$ =  0.83392334

$q_3$ =  0.07568359

The response follows the ideal -3 dB/octave curve to within $\pm$0.3 dB error over a 10 octave range from 0.0009$\times \pi$ to 0.9$\times \pi$.  ("Nyquist" = $\pi$.)


Principal Component Analysis (PCA) on Convolutional Neural Network (CNN) Features Please, I have a question regarding PCA and features which are extracted from a convolutional layer. link

if we have a test dataset , and we extract all conv features of all images at test dataset called feat. Then we do the following normalization and PCA to all conv features

feats = normalize(feats)
pca = PCA(512 ,whiten=True)
pca.fit(feats)

And then use pca model for test dataset

# PCA MODEL - use paris for oxford data and vice versa
self.pca = pickle.load(open('_oxford.pkl', 'rb'))

And then applying PCA on the test dataset using the previous pca

print "Applying PCA"
self.db_feats = normalize(self.db_feats)
if self.pooling is 'sum':
self.db_feats = self.pca.transform(self.db_feats)
self.db_feats = normalize(self.db_feats)

What is the idea of using pca model of another dataset to transform the features of test dataset?

I'm not into details of this specific case but I can see some logic.

A convolution layer can be reformulated as a Matrix Multiplication:

$$ y = W x $$

Let's say we trained on Data Set $ {x}^{1} $ which is big and general.
Namely we expect the trained weights $ {W}^{1} $ to be good enough for almost any other data set.

Let's assume we have another data set $ {x}_{2} $ with its matching weights $ {W}^{2} $.
By our assumption about how good $ {x}^{1} $ generalizes it means $ {W}^{1} \sim {W}^{2} $ with regard to some metric.

If this metric is the Frobenious Norm / Matrix Induced $ {L}_{2} $ Norm then using PCA (Which is basically SVD) makes sense as using the PCA we take the component of $ {W}^{1} $ which generalizes the most and use it for other case.


Amplitude vs Frequency Response Consider the simple first-order low pass filter, described by

$\tau \frac{d \mathcal{O}}{dt} + \mathcal{O}(t)=\mathcal{I}(t)$

Considered as a linear system (ie $\mathcal{L} \mathcal{O} = \mathcal{I}$), $\mathcal{L}$ has transfer function $\frac{1}{1 + i \omega \tau}$ (up to maybe some constants, minus signs etc etc depending on convention). By transfer function, I mean $\tilde{R}(k)$ such that if $\tilde{\mathcal{I}}$ denotes the Fourier/Laplace transform of the input, then $\tilde{\mathcal{O}}(t) = \tilde{R}(t)\tilde{\mathcal{I}}(t)$. Questions:

1) Is the transfer function the 'same thing' as the frequency response? The reference textbook I have (Signals and Systems by Oppenheim & Willsky) talks of the frequency response of the 1st order LPF $H(i\omega) = \frac{1}{1 +i \omega \tau}$ which as far as I can tell amounts to something very similar.

2) If we solve the LPF ODE for a sinusoidal input $sin(\omega t)$, then we get an expression for the output amplitude $A(\omega) = \frac{1}{\sqrt{1+\omega^2 \tau^2}}$. How is this related to the transfer function/frequency response?

Supplementary information/question: I'm actually a mathematics student in the situation of having to learn a lot of signals processing very quickly this summer. Are there any textbooks written out there which explain signals and systems in the perspective of a maths student? As far as I can tell, most books seem to pitched towards engineers (and also use $j$ instead of $i$), and are a little different to what I'm used to - I don't want to learn lots of things twice from two different perspectives.

Thank you very much for you time.

The answer to your suplemental question: watch Openheim's free youtube digital signal processing MIT lectures, especially the lectures number 2,3,4, and you will get answers on the first question as well. The lectures are well grounded in math, and I am sure that this is the most efficient way to learn DSP on your own, from the point of view of a mathematician.


Noise rejection / feature extraction I'm hoping someone could point me in the direction of some signal processing methods to clean up my data. I'm collecting physiological data from human muscle (lower leg; gastrocnemius) during walking/gait.

I know the sensor is picking up contractions, but the problem is it's so plagued by motion and skin artefacts I can barely do anything with it. Performing an isometric contraction (contracting the muscle without changing its length) shows a good signal (as there are no motion or skin artefacts to contaminate it. However, when collecting real-world data of a subject walking I'm getting huge spikes on the heel strike, amongst other unpleasantries.

The sensor type revolves around a MEMS microphone, so basically it picks up everything. It is un-amplified and only filtered by post processing. The sampling rate is 1kHz and seems to be more than enough. It's a single channel/vector of a continuous signal (i.e. 10 minutes of collection at 1kHz is a 1x600000 vector).

Please help me by suggesting some methods on what I can do. I've been looking at wavelet transforms which look promising, but any help on the matter will be much appreciated. Also look at PCA but having a vector makes it difficult (unless I break my data into windows?). I'm trying to keep it as uncontrolled as possible, so any analysis that I can do on my current collected data would be great!

Can provide some data if it helps in your analysis.

Cheers!

I've tried a copule of methods on your data (low-pass filtering via DFT, SSA filtering), but I don't know what result will satisfy you.

You mentioned that the signal is sampled with a higher-then-enough frequency. Can you rise the  sampling frequency even higher? I'm talking about oversampling. This can help to get rid of some part (some kind) of the noise present.

You also mentioned PCA and that it's not appropriate for univariate time series. Right, PCA is almost useles with univariate time series, but there is an appropriate and closely related technique, SSA. You can find detailed description here: Singular Spectrum Analysis for Time Series (free sample is enough: Chapter 2: Basic SSA). And my own GNU Octave (hope Matlab compatible) implementation here: ssa-octave.m (there is also Scilab implementation, but you probably don't need it). Take a look at chapter 2.4 "Choice of Parameters in Basic SSA" for interpretation and tuning tips.

It's also not clear to me: the "5 seconds of stationary" pattern in isometric data sample should be filtred out or it is a useful signal too?


Any good C# library for signal processing? Does anyone can link me to a good C# DSP library? I saw some good C++ libraries but I don't have much with C++. Is there any way to wrap them ?

My main porpuse is to apply some filters(Lowpass, Highpass,Notch....) on a a buffer. Thanks

I can recommend these libraries:

BASS - also has lot's of functions. Although it is not written for C#, they provide you with a wrapper (Bass.Net). Keep in mind that if you want to use this library in commercial project then you have to pay for the license.

NAudio - is a mature project with lots of functions and support for variety of codecs.

cscore - rather new project with some nice capabilities, especially for visualizations, and some codecs that NAudio is missing. You might want to give it a try at some point, or even contribute to it.


Aliasing in wave-scanning oscillator I wrote a wavetable oscillator that scans a file and extracts grains etc. The sound doesnt seem clear and I think it has aliasing.

What are ways to be able to play the audio file at any frequency without aliasing?

I know there should be interpolation, is there any specific type of interpolation that is better suited for recorded audio files?

Are there any other techniques to reduce this aliasing during playback?

I am unable to visualize how an audio file at any frequency, that is played back at different frequencies by jumping samples during indexing, aliases back into the waveform, and what depth of interpolation is needed.

Will scanning an audio file at higher frequencies cause aliasing due to index jumping, or is it rather precision loss? If the audio is already recorded, without aliasing, 44khz sample rate, bandwidth of 22khz and playback rate up to +6 octaves, is the aliasing added inside the real time signal?

i am a little confused by the first sentence. i dunno how a wavetable oscillator scans a sound file at all. some algorithm can scan a file and extract either grains (or little "wavelets", which is a term i have misused in the past) or single-cycle wavetable "snapshots". but the wavetable oscillator starts with wavetable(s) and outputs a periodic or quasiperiodic tone.

about aliasing, if you extract wavetables from a real tone (and let's assume you have a decent pitch detector, so you know how long each period is to a precision of fraction of sample, and decent interpolation, so you you can resample the period to a, say, 256 or 1024 point wavetable) and you up-pitch the tone by more than a couple of semitones, there is a danger of aliasing unless you remove some of the higher harmonics. you can remove these higher harmonics by performing a DFT (a.k.a. FFT) on the wavetable and zero some of the bins centered around bin N/2 . then inverse DFT back to a wavetable.

actually, the number of bins centered around bin N/2 that are already zero is dependent on the degree of upsampling in the interpolation when the period in the soundfile is converted into a wavetable.  consider this example:  suppose the sample rate of the sound file is 44100 Hz and you have a note that is at middle C, with a fundamental at 261.6255653 Hz.  this means that your period length is 168.5615 samples and a really good pitch detector can return a period value with that precision (even though i doubt it would be that accurate, precision and accuracy are different things).  now a 128-sample wavetable would be insufficient to represent that waveform, or at least possibly insufficient to represent all of the harmonics of that waveform.  there are conceptually 84 harmonics.  a 256-samplewavetable would be sufficient to represent all of the harmonics, but let's say, for playback reasons, you resample that 168.5615 sample period to a 1024-sample wavetable, $x[n]$.  that means you have a step size (the reciprocal of the sample-rate-conversion ratio) of 168.5615/1024 = 0.16461084 .

now, if you were to FFT that 1024-sample wavetable, you would possibly have non-zero bins for $X[0]$ up to $X[84]$ (and also the mirror image of $X[1024-84]$ up to $X[1023]$ .  all of the bins from $X[85]$ to $X[1024-85]$ should be zero if the interpolation kernel was any good.

now, when you play that wavetable back at the same note, middle C, then the sample stride would be the reciprocal of the step size above, or 6.074934 wavetable samples per output sample (assuming the same 44100 Hz sample rate).  and you would have no aliasing because all of those harmonics above the 84th harmonic are zero.  but if you start pitching that note higher, your stride will be longer.  if you pitched it 7 semitones higher (to G), then the stride would be 9.1021167 wavetable samples per output sample.  but then your 84th harmonic gets bumped up to around 33 kHz which is above Nyquist and you will have aliasing.

so what you need to do, if you plan to pitch that note up 7 semitones (i would not pitch a sampled piano up that far, because it will sound "plinky"), you must lose all of those harmonics above the 56th harmonic.  ($56 \approx 84 \times 2^{-7/12} $).  you do that with the FFT, zeroing bins, and then inverse FFT