Distorted output of delay VST effect in C I'm getting some odd results with some code I'm working on. Help me with the implementation of this delay effect?

The code takes an input, adds delayed input to it and then passes it to the output whilst also copying it back to the delay line. Control of the delay level and feedback are provided to control respectively the amount of delayed signal added to the output and returning to the delay line.

Here's the process of my delay effect.

void delay::processReplacing(float** inputs, float** outputs, VstInt32 sampleFrames)
{
        float *in = inputs[0];
        float *out = outputs[0];
        float delayed_sample;
        int delayed_sample_index;
        for (int frame_idx=0;frame_idx< sampleFrames;frame_idx++)
        {
                delayed_sample_index = (((int)(cBuf.index-delay_samples) + CIRCULAR_BUFFER_LENGTH) % (int)CIRCULAR_BUFFER_LENGTH);
                delayed_sample = (cBuf.buffer)[delayed_sample_index];
                //write current sample to buffer
                (cBuf.buffer)[cBuf.index] = in[frame_idx] + (delayed_sample*feedback);
                //Calculate + write output sample
                out[frame_idx] = in[frame_idx] + (delay_level*delayed_sample);
                //Increment buffer index
                cBuf.index = ((cBuf.index+1)) % CIRCULAR_BUFFER_LENGTH;
        }
}

On the whole it works pretty well and functions as a delay effect should. However, the output I'm getting is slightly distorted. I'm not sure if this is an error in the arithmetic of sample times and values or something else. Any ideas?

I am not familiar with signal processing for audio, but I would like give answer which contain common programming suggestions. Sorry, I really do not know your programming level and may be this suggestions will be look obvious and ridiculous, but I hope it will be useful.

  I'm getting some odd results with some code I'm working on.

When you get odd results usual this caused by one of three reasons:

incorrect implementation of correct logic;
correct implementation of incorrect logic;
incorrect implementation of incorrect logic.

Some implementation errors such as: using of uninitialized variables, using of wrong pointers, accessing to memory out of bounds can be easy found when using some special tools. One kind of tools is static code analyser, there are a large number of varied of static code analysers, but I would suggest you use cppcheck. Other kind of tools is (memory) profilers. I would suggest you use memcheck tool from valgrind suite. Very often this tools helps find very tricky errors.

To find logical error you need test code with different input data and conditionals. It is possible to do it manually time to time, but more easy to do it automatically by calling small code snippets. This code snippets can be designed as simple console application or as test case using some testing framework. Testing framework is not necessary for writing test, but it need for manage tests. I suggest to use the C++ Google testing framework.

Now I will try to implement my own suggestions.

First of all to allow a piece of existed code will be compiled without errors I added header file and put into it all required types and variables:

#include <cstdint>

using VstInt32 = int32_t;
constexpr int delay_samples             {1};
constexpr int CIRCULAR_BUFFER_LENGTH    {4};
constexpr float feedback             {0.0f};
constexpr float delay_level          {0.0f};

struct CircularBuffer
{
    CircularBuffer ()
    {
        for (int i = 0; i < CIRCULAR_BUFFER_LENGTH; ++i) {
            buffer [i] = 0.0f;
        }
        index = 0;
    };

    int index;
    float buffer [CIRCULAR_BUFFER_LENGTH];
};

class delay
{
    public:
        void processReplacing (float **, float **, VstInt32);

    private:
        CircularBuffer cBuf;
};

Shared code I placed into source (*.cpp) file:

#include "delay.h"

void delay::processReplacing (float ** inputs, float ** outputs, VstInt32 sampleFrames)
{
    float * in  = inputs  [0];
    float * out = outputs [0];
    float delayed_sample = 0.0f;
    int delayed_sample_index;
    for (int frame_idx = 0; frame_idx < sampleFrames; frame_idx++) {
        delayed_sample_index = (((int)(cBuf.index - delay_samples) + CIRCULAR_BUFFER_LENGTH) % (int)CIRCULAR_BUFFER_LENGTH);
        delayed_sample = cBuf.buffer [delayed_sample_index];
        // Write current sample to buffer
        cBuf.buffer [cBuf.index] = in [frame_idx] + (delayed_sample * feedback);
        // Calculate + write output sample
        out [frame_idx]          = in [frame_idx] + (delayed_sample * delay_level);
        // Increment buffer index
        cBuf.index = (cBuf.index + 1) % CIRCULAR_BUFFER_LENGTH;
    }
}

Now time to write tests. I am not familiar with logic of function which you try to implement, so I write one simple test just for example:

#include <gtest/gtest.h>
#include <iostream>
#include "delay.h"

class delay_test : public ::testing::Test
{
    public:
        void SetUp    () override;
        void TearDown () override;

        static constexpr int sampleFrames {8};
        float ** inputs;
        float ** outputs;
};

// This code will be automatically called before every test.
// This is feature of Google Tests Framework.
void delay_test::SetUp ()
{
    inputs      = new float * [1];
    outputs     = new float * [1];
    inputs  [0] = new float   [sampleFrames];
    outputs [0] = new float   [sampleFrames];
    for (int i = 0; i < sampleFrames; ++i) {
        inputs  [0][i] = 0.0f;
        outputs [0][i] = 0.0f;
    }
}

// This code will be automatically called before every test.
// This is feature of Google Tests Framework.
void delay_test::TearDown ()
{
    delete [] inputs  [0];
    delete [] outputs [0];
    delete [] inputs;
    delete [] outputs;
}

TEST_F (delay_test, response_on_1_Success)
{
    // Arrange.
    delay d;
    inputs [0][0] = 1.0f;

    // Act.
    d.processReplacing (inputs, outputs, sampleFrames);

    // Assert.
    for (int i = 0; i < sampleFrames; ++i) {
        if (1 == i) {
            EXPECT_FLOAT_EQ (1.0f, outputs [0][i]);
        }
        else {
            EXPECT_FLOAT_EQ (0.0f, outputs [0][i]);
        }
    }
}

Last step is writing build script. I use cmake, so simple build script has the following content:

cmake_minimum_required (VERSION 2.8)

project (delay)

set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -mtune=native -std=c++11 -Wall -Wextra")

set (
    delay_test_SOURCES
    ${PROJECT_SOURCE_DIR}/delay.cpp
        ${PROJECT_SOURCE_DIR}/delay_test.cpp
)

include_directories (
    ${PROJECT_SOURCE_DIR}
)

add_executable (
    delay_test
    ${delay_test_SOURCES}
)

target_link_libraries (
    delay_test
    gtest
    gtest_main
    pthread
)

Now it is possible to test code.

Let's check code with cppcheck:

$ cppcheck -I. --enable=all *.cpp
Checking delay.cpp...
[delay.cpp:7]: (style) The scope of the variable 'delayed_sample' can be     reduced.
[delay.cpp:8]: (style) The scope of the variable 'delayed_sample_index' can be reduced.
1/2 files checked 33% done
Checking delay_test.cpp...
2/2 files checked 100% done
Checking usage of global functions..
[delay_test.cpp:19]: (style) The function 'SetUp' is never used.
[delay_test.cpp:33]: (style) The function 'TearDown' is never used.
(information) Cppcheck cannot find all the include files (use --check-config for details)

The results are good, nothing suspicious is found.

Now let's launch tests and simultaneously memory profiler:

$ valgrind ./delay_test
==11407== Memcheck, a memory error detector
==11407== Copyright (C) 2002-2013, and GNU GPL'd, by Julian Seward et al.
==11407== Using Valgrind-3.10.0 and LibVEX; rerun with -h for copyright info
==11407== Command: ./delay_test
==11407==
Running main() from gtest_main.cc
[==========] Running 1 test from 1 test case.
[----------] Global test environment set-up.
[----------] 1 test from delay_test
[ RUN      ] delay_test.response_on_1_Success
/home/gluttton/circ_buf/delay_test.cpp:57: Failure
Value of: outputs [0][i]
  Actual: 1
Expected: 0.0f
Which is: 0
/home/gluttton/circ_buf/delay_test.cpp:54: Failure
Value of: outputs [0][i]
Actual: 0
Expected: 1.0f
Which is: 1
[  FAILED  ] delay_test.response_on_1_Success (81 ms)
[----------] 1 test from delay_test (95 ms total)

[----------] Global test environment tear-down
[==========] 1 test from 1 test case ran. (144 ms total)
[  PASSED  ] 0 tests.
[  FAILED  ] 1 test, listed below:
[  FAILED  ] delay_test.response_on_1_Success

 1 FAILED TEST
==11407==
==11407== HEAP SUMMARY:
==11407==     in use at exit: 0 bytes in 0 blocks
==11407==   total heap usage: 306 allocs, 306 frees, 61,494 bytes allocated
==11407==
==11407== All heap blocks were freed -- no leaks are possible
==11407==
==11407== For counts of detected and suppressed errors, rerun with: -v
==11407== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)

Good news is that no errors was found by valgrind (we can realize it by valgrind output messages). But bad news is that test is fell, it means that code works not as I expect.

Conclusions: in my opinion your code implemented without visible implementation errors, but you need check logic of your code with various input data.


How to Linearly Combine Two Unbiased Estimators of One Parameter without Knowledge of Their CoVariance? I have two unbiased estimators of one parameter, $\tau$. The first estimator, $r_1$, is the better estimator with lower variance than the second estimator, $r_2$.

I also have:
$ \mathbb{E} \left[ {r}_{1} \right] = 0.8 \tau $ , and $ \mathbb{E} \left[ {r}_{2} \right] = 0.2 \tau $.

I want to use both estimators to achieve lower variance. Is there a systematic way to find a linear combination of $r_1$ and $r_2$ that results in a lower variance and unbiased estimate, without a knowledge of the variance and covariance of the estimators?

Basically I'd do something like:

$$ \frac{ {\sigma}_{2}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2}} {r}_{1} + \frac{ {\sigma}_{1}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2}} {r}_{2} $$

This is the optimal weighing given knowledge of the Variance only (Well, linear).
Basically it assumes the cross correlation is 0.

Derivation

There are many ways to derive this.
The easiest would be intuition of knowing that Variance is the Linear Measure of the quality of an estimator. But we're leaving intuition our for this moment.

So easiest formal derivation could be using the ideas @Mark Leeds talked about.

Our model is:

$$ {r}_{3} = {w}_{1} {r}_{1} + {w}_{2} {r}_{2} $$

We want to guarantee the mean of $ {r}_{3} $ is un biased so:

$$ \mathbb{E} \left[ {r}_{3} \right] = \mathbb{E} \left[ {r}_{1} \right] = \mathbb{E} \left[ {r}_{2} \right] = {w}_{1} + \mathbb{E} \left[ {r}_{1} \right] + {w}_{2} + \mathbb{E} \left[ {r}_{2} \right] \Rightarrow {w}_{1} + {w}_{2} = 1 $$

So we'll use $ {w}_{2} = 1 - {w}_{1} $. Now we want to minimize its variance. Since it is not biased we're after minimizing its second moment:

$$\begin{align*}
\mathbb{E} \left[ {r}_{3}^{2} \right] & = \mathbb{E} \left[ { \left( {w}_{1} {r}_{1} + \left( 1 - {w}_{1} \right) {r}_{2} \right) }^{2} \right] \\
& = {w}_{1}^{2} \mathbb{E} \left[ {r}_{1}^{2} \right] + { \left( 1 - {w}_{1} \right) }^{2} \mathbb{E} \left[ {r}_{2}^{2} \right] + 2 {w}_{1} { \left( 1 - {w}_{1} \right) } \mathbb{E} \left[ {r}_{1} {r}_{2} \right] \\
& = {w}_{1}^{2} {\sigma}_{1}^{2} + {\left( 1- {w}_{1} \right)}^{2} {\sigma}_{2}^{2}
\end{align*}$$

We assumed $ \mathbb{E} \left[ {r}_{1} {r}_{2} \right] = 0 $.

Now, this is a quadratic form so to minimize just derive it and find what vanishes it:

$$ {w}_{1} {\sigma}_{1}^{2} - {\sigma}_{2}^{2} + {w}_{1} {\sigma}_{2}^{2} = 0 \Rightarrow {w}_{1} = \frac{ {\sigma}_{2}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2} }, \; {w}_{2} = \frac{ {\sigma}_{1}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2} } $$

Pay attention that the solution doesn't require us to know the exact variances of the estimators being fused. Only their ratio.

So in the case above where $ {\sigma}_{1}^{2} : {\sigma}_{2}^{2} = 0.8 : 0.2 $ the solution is $ {w}_{1} = 0.2 $ and $ {w}_{2} = 0.8 $. This is the best you can do without the knowledge of the co variance.

This fused veritable will have variance $ {\sigma}_{3}^{2} = {0.2}^{2} {\sigma}_{1}^{2} + {0.8}^{2} {\sigma}_{2}^{2} = 0.16 \hat{r} $

Wikipedia Solution

You may, as a reference, look at Wikipedia - Sensor Fusion.
They wrote the same equation in a real clumsy manner:

Yet they are equivalent (Just multiply their result by $ \frac{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2} }{ {\sigma}_{1}^{2} + {\sigma}_{2}^{2}} $)

Pay attention that the sensor model assumes that 2 sensors are not correlated (Well, it assumes noise is property of the measurement equipment, hence this assumption is made).


Can white Gaussian noise improve signal level in non-linear device? In the system shown below, can white Gaussian noise amplify or attenuate the input signal?

Most of the time we check signal level proportional to noise level (such as noise figure) not singly.

I have obtained PSD of my desire system's output and it is a combination of

a) Dc term
b) reaction of input noise to itself
c) reaction of input signal to itself
d) reaction of input signal and input noise.

At system's output at frequency $\omega$ should I assume terms (c) and (d) as the signal level? Or c is enough?

As Jason R. says in the comment, dither will do this.  Another phrase to look for is noise shaping where the noise spectrum can be reshaped to not interfere with the signal of interest.


IIR filter parallelization Is there any current ideas about IIR filter algorithm parallelization with thread count greater than filter order?

Ideally, parallelization method should load GPU efficiently.

You might be interested in the following paper:

Nehab, D; Maximo, A; Lima, R; Hoppe, H: GPU-Efficient Recursive Filtering and Summed-Area Tables. ACM Trans. Graph. 30(6):176 (December 2011). doi://10.1145/2024156.2024210
http://w3.impa.br/~diego/publications/NehEtAl11.pdf

The basic idea is to recognize that linear filters are associative (at least in theory) and thus you can

divide your input into chunks.
every thread initializes its iir filter state with 0s and processes its chunk.
every thread stores the final state of its iir filter for its chunk.
now (because of associativity) you can figure out what the beginning filter state is of every chunk.
so then you run the filter over the chunk again, this time correctly initialized and get the correct output.

As @jonsca pointed out, you were a bit vague about your use case.  The above all assumes that you can do batch processing on your input.  If you need some kind of real-time filtering none of what I said is relevant.


Applying FFT to PCM stream I would like to use a .WAV file and convert it into a spectogram. Any idea how to do this with Naudio? Basically the flow should be like this, I think:

Convert WAV file to PCM stream
Convert PCM stream to double-array and applying FFT
Using the result from (2) and convert it to a spectrogram

Your description of how to create a spectrogram is correct. Here's an explanation with a little more detail:

First, create a WaveChannel32, using a WaveFileReader. This will output a PCM byte stream. Then read 4 times the size of your FFT bytes from the WaveChannel32 into a buffer. Convert every four bits in the buffer into an integer.

It's important at this point to use a window function on our interval, to emphasise amplitudes closer to the current time stamp.

I'd then suggest using the Exocortex FFT library to calculate the FFT. The Fourier.FFT method will convert an array of complex numbers to its corresponding Fourier transform. Then ignore the second half of the array produced, to prevent aliasing.

The first half of the array corresponds to the column of the spectrogram at our current time stamp. Each number in the array corresponds to a single pixel in the spectrogram, the frequency of which is determined by the index of the array (equally spaced throughout 0 to 22050Hz). The darkness of each pixel is determined by the magnitude squared of the complex number corresponding to said pixel.

Here's some sample code for producing one column:

var waveChannel = new NAudio.Wave.WaveChannel32(new NAudio.Wave.WaveFileReader("file path"));

var buffer = new byte[4096*4];

waveChannel.Read(buffer, 0, 4096*4);

var samples = new ComplexF[4096];

for (var i = 0; i < 4096, i++)
{
    samples[i].Re = BitConverter.ToSingle(buffer,i*4);

    // Window should be pre-populated with values from your window function
    samples[i].Re *= Window[i];
}

Fourier.FFT(samples, FourierDirection.Forward);


Finding a signal inside another signal I have a 1D array of 32 bit floating point numbers representing the current draw of a radio operating over a 40 minute span.  The number of elements in the array is 29 million and change.  There is a signal of interest inside this array that I would like to be able to identify in future signals.  I've copied the part of the array that I am interested in and used a python script to generate a parameterized mask.  I would like to use this mask to look for similar signals in my 40 minute array.  The issue is the amount of time it takes to "sweep" this mask over all the data to see if all the data points fit under the mask.

Generally, I would like to be able to identify signals of interest in large data sets and extract their indices.

I am using Python, numpy, and SciPy.

Images 1 Small signal and mask 2 Full 40 minute radio current draw

Given a signal waveform of interest, what is the most effective way to detect whether that signal exists inside another signal?

Here is a link to the raw data.  It's a zipped csv format.  The values are time and current draw in amps.  Sampling rate of the data:  One sample every 8.192E-05 seconds.  Signal of interest begins around 917.65 seconds and lasts less than one second.

https://www.dropbox.com/s/rw5yt5zv5rpraoh/ExportedData.zip?dl=0

As I wrote in the comments, this can be solved using cross-correlation.

Take your mask of length $L$, invert it and start performing the convolution between this inverted mask and the original signal (or, directly, use scipy.signal.correlate() in Python with the two signals, the original one and the mask). If you get a peak in your correlation at sample $m$, this means that if you take your mask and overlay it over the original signal such that the first sample of the mask is at the $m$-th sample of the original signal (using the function with the parameter mode='same'), then you will have a good match (i.e., for that instant of time, the signal is pretty similar to the mask). Then you can conclude that from $m$ to $m+L−1$ the signal has a shape that looks like the one of the mask.

I checked the signal you've uploaded and, as you pointed out, following the procedure above one gets some non-desired peaks due to the fact that the original signal has some high values that don't correspond to a similarity in shape with the mask, but nevertheless return a high value of cross-correlation because the sum of the elementwise product is high. This could be solved "filtering out" these non-desired values.

The most similar signal to the mask will be... the mask. Thus, calculate the correlation between the mask and itself (the autocorrelation at $0$). This will give you a number that would represent the maximum possible similarity to the mask. Then, you can ignore those values in the original cross-correlation that are greater than the autocorrelation at $0$, as they surely correspond to those high peaks the original signal has. This should do the trick.


Common symbols or standards used for DSP block diagrams? When creating high level block diagrams of DSP or signal processing, what symbols are typically used?
Are there any standards? For example its common to see gain to be a triangle. Whilst for a low pass filter, I have seen a few ways.
Assuming an absence of standards, whats a practical way of finding/referencing symbols? Was thinking there must be some reasonable comprehensive software packages, say like Simulink, that would have made an intelligent attempt at creating symbols. Or worse, review various DSP books? Can anybody recommend an approach or resource? In an ideal world, an open source project set of .svgs!

(source: researchgate.net)

The best option would be to review the standards published by large standards organisations such as ANSI, ISO, DIN, EN and others. Other countries have respectable standards organisations but these are the ones I happened to have off the top of my head at the moment.

Usually however, you might find that there is a price associated with getting access to a particular standard. Having said that...

  When creating high level block diagrams of DSP or signal processing, what symbols are typically used?

I think that what has the potential to cover everything is the flowgraph, which is not too far off the underlying abstract tool for representing interconnected components.

More specific to DSP are only a handful of symbols. These are, the delay element (a square with $z^p$ where $p$ denotes the delay and its sign is usually negative), the "summation node" (a circle with a $+$ or $\Sigma$ inside it) and the "gain" which denotes multiplication by some coefficient (so, really, a circle with some number inside it).

I have quoted the "summation node" above because you might find diagrams where you have another operation taking place at a similar node. For example, you might have "modulation" (the circle with the X cross occupying the full diameter of the circle).

  Are there any standards? For example its common to see gain to be a triangle.

The triangle is a reference to the symbol for the amplifier and I think that it makes obvious sense to everyone who has studied electronics.

This is another point: Diagrams have a context. On a CPU, the term "coefficient" or "multiplication" makes more sense than "amplification" or "gain". If you are talking to Radio Frequency (RF) people, "multiplication" and "amplification" are more likely to be associated with "mixer" and "amplifier" in the real world, but if you are talking to an Audio person and you tell them about "mixing", they might think something quite different.

  Whilst for a low pass filter, I have seen a few ways.

As far as I can tell from an old stencil ruler, the wavy representation of filters is a DIN standard and you strike out the cut off bands. So, a band pass filter has three waves with the low and high crossed out.

I have doubts for that demodulator though because it denotes an envelope detector which means AM demodulation. If you left "demodulator" open, without definition, I think someone might associate it with a "mixer".

Hope this helps.


how can I design Chebyshev filter with 5th-order polynomial and 0.001%-allowed ripples? I'm very beginner in signal processing. I need to design a Chebyshev filter with 5th-order polynomial and 0.001%-allowed ripples by using matlab as I found in a paper. I'm filtering pressure data inside combustion engine. In matlab I've to define normalized passband edge frequency and peak-to-peak passband ripple.
I don't know what did the paper's author mean by 0.001% allowed ripples?

Any help would be appreciated

A type I Chebyshev filter has a passband magnitude response oscillating between the values $1$ and $1-\delta$, where $\delta$ is the maximum passband approximation error. The Matlab function $\tt{cheby1}$ requires the desired passband ripple in dB:

$$R_p=20\log\left(\frac{1}{1-\delta}\right)$$

The normalized cut-off frequency is

$$\omega_c=2\pi\frac{f_c}{f_s}$$

where $f_c$ is the cut-off frequency in Hertz, and $f_s$ is the sampling frequency in Hertz. The input parameter $\tt{Wc}$ required by $\tt{cheby1}$ is equal to $\omega_c/\pi$.

I don't think that the authors of the paper you refer to got the percentages right. A tip: for a "ripple percentage" of $0.001$ try $\delta=0.0001$.


Volterra Kernel Convolution Method Here I understand that the first term is a simple convolution of the input signal with first volterra kernel that I have acquired through farina sweep method, my question is if the second term is the convolution of second volterra kernel with our input signal raised to the power of two? Would the code using matlab, of the first and second terms be correctly written as `

outputsignal1=conv(voltkernel1,inputsignal.^1); %basic convolution with system impulse
outputsignal2=conv(voltkernel2,inputsignal.^2); %convolution of second kernel with squared signal
summedsignal=outputsignal1+outputsignal2;  %dynamic nonlinear output signal with memory

Yes the MATLAB code is correct. Be careful though, the bandwidth of the signal squared is twice that of the signal itself, which may lead to aliasing if the sampling frequency is too low compared to the signal bandwidth. This can be remedied by properly resampling the signal to a higher sampling frequency with a lowpass filter, then squaring, low-pass filtering, and decimating back to the original rate.

Another pitfall if you are doing everything by yourself is that the impulse-response-like partitions of the deconvolved sweep response are not directly the Volterra kernels. The rest of this answer will discuss how to get usable kernels from the partitions corresponding to the harmonics.

First, you may need to phase shift each partition with a constant phase shift, depending on what kind of sweep you used. You can test if there is any phase shift needed (and test if you do the phase shift correctly), by using as a test system a memoryless non-linearity that is the sum of a number of early Chebyshev polynomials of the first kind, as this will generate for a unit amplitude cosine input a number of harmonics in the same phase (explanation further down in this answer). If I understood the following paper correctly, for certain kinds of sweeps there is no need for phase shifts:

Antonín Novák, Laurent Simon, František Kadlec, Pierrick Lotton. Nonlinear System Identification Using Exponential Swept-Sine Signal. IEEE Transactions on Instrumentation and Measurement Volume: 59, Issue: 8, Aug. 2010.

Then each one-dimensional Volterra kernel can be calculated as the weighted sum of the (phase-shifted) partitions, with the weights obtained from solving a simple linear system. This is covered in Angelo Farina's papers, but at least some of them complicate the derivation by using sines rather than cosines, which gets into the earlier phase shift territory. Cosines are simpler because they have real frequency-domain representations implying no phase shift. Here is an Octave script (using cosines) that solves the weights given the number of harmonics:

N = 3; #Number of harmonics including the fundamental (you can touch this)
M = 2*N+1; #Matrix height, must fit positive, zero and negative frequencies
x = zeros(M, N+1);
# From cosines to complex exponentials. Negative indexing by shifted columns
x(1, 1) = 1; #DC bias
x(3, 2) = 0.5; #This is negative fundamental
x(1, 2) = 0.5; #This is positive fundamental
for k = 3:N+1
  x(:, k) = conv(x(:, k-1), x(:, 2))(1:M); #Time domain exponent k
endfor
for k = 1:N+1
  x(:, k) = shift(x(:, k), 1-k); #Shift back columns
endfor
x = x(1:N+1, :); #Go back to cosines
x(2:N+1,:) .*= 2 #
inv(x) #Solve system

The output is:

x =
   1.00000   0.00000   0.50000   0.00000
   0.00000   1.00000   0.00000   0.75000
   0.00000   0.00000   0.50000   0.00000
   0.00000   0.00000   0.00000   0.25000

ans =
   1  -0  -1  -0
   0   1  -0  -3
   0   0   2  -0
   0   0   0   4

The answer ans is the inverse of matrix x. The answer means that if you have three harmonics including the fundamental, then you go from the partitions $h_n$ corresponding to harmonics $n \in \{1,\, 2,\, 3\}$ to kernels $k_m$ corresponding to exponents $m \in \{1,\, 2,\, 3\}$ by:

$$\left\{\begin{eqnarray}k_0 &=& h_0-h_2\\
k_1 &=& h_1 -3h_3\\
k_2 &=& 2 h_2\\
k_3 &=& 4 h_3\end{eqnarray}\right.$$

This is the solution of:

$$\begin{align}&k_0 + k_1 \cos(t) + k_2 \cos^2(t) + k_3 \cos^3(t) = h_0 + h_1 \cos(t) + h_2 \cos(2t) + h_3 \cos (3t)\\
\Rightarrow & k_0 + \frac{k_2}{2} - h_0 + \left(k_1 + \frac{3k_3}{4} -h_1\right)\cos(t) + \left(\frac{k_2}{2} - h_2\right)\cos(2t) + \left(\frac{k_3}{4} - h_3\right)\cos(3t) = 0,\end{align}$$

where $h_0$ is a zero frequency bias that is not determined by the sweep method and is normally assumed zero, $h_0 = 0,$ and $k_0$ is a term that eliminates the bias. I'm not sure how to use it, probably just add to each output sample the sum of samples of $k_0.$

The solution matrices for smaller N are contained in the solution matrices for larger N. The matrices follow the On-line Encyclopedia of Integer Sequences (OEIS) sequence A053120. OEIS gives simpler ways to calculate the weights, and also hints to a link to Chebyshev polynomials of the first kind. They have coefficients that are identical to the weights. The link is actually such that if you use as the memoryless non-linear elements not powers but Chebyshev polynomials of the first kind, $T_n(x)$, then the kernels to use are directly the (phase-shifted) partials corresponding to the harmonics. This comes from the property $T_n\big(\cos(t)\big)=\cos(nt).$ This is also covered in:

Antonin Novak, Laurent Simon, Pierrick Lotton, Joel Gilbert. Chebyshev Model and Synchronized Swept Sine Method in Nonlinear Audio Effect Modeling. Proc. of the 13th Int. Conference on Digital Audio Effects (DAFx-10), Graz, Austria, September 6-10, 2010.

The N = 7 solution from the script agrees with:

Marc R´ebillat, Romain Hennequin, Etienne Corteel, Brian Katz. Prediction of harmonic distortion generated by electro-dynamic loudspeakers using cascade of Hammerstein models. 128th Convention of the Audio Engineering Society, May 2010, London, United Kingdom. pp.7993, 2010. <hal-00619357>


Verifying the computation of a convolution I have an input signal $$x(n)=\left(3,-5,4,3,-1,-2,6,8\right), n=-3,..,4$$ and impulse response $$h(n)=(1,-1,1,-1,1), n=-1,...,3.$$

The convolution between $x(n)$ and $h(n)$ is

$$x(n)*h(n)=\sum_{-\infty}^\infty x(k)h(n-k)$$
If I'm not mistaken, I can reduce this to the finite sum
$$=x(-3)h(n+3)+x(-2)h(n+2)+x(-1)h(n+1)+x(0)h(n)+x(1)h(n-1)+x(2)h(n-2)+x(3)h(n-3)+x(4)h(n-4)$$

So, if I were to calculate the sum by hand, I would proceed by first evaluating

$x(-3)h(n+3)=3h(n+3) = (3,-3,3,-3,3)$ for $n=-4,...,0$
$x(-2)h(n+2)=-5h(n+3) = (-5,+5,-5,+5,-5)$ for $n=-3,...,1$
etc...

Then I would only add up the terms whose position $n$ aligns (that is, add all terms at n=-4 together, all terms at n=-3 together, etc...).

Is my approach to evaluating the convolution correct?

I would have liked to give a longer answer, but yes, your approach is correct. I don't see any problems with what you laid out in your question.


Why do we need to estimate eigenvalues? I am not working in signal processing field, but recently I happen to read a paper which estimates source numbers using Gerschgorin radii, and I feel kind of confused about why we need to estimate eigenvalues (this might be a silly question in your point of view). I google it and some answers say the reason is that it is hard to get eigenvalues in some cases. But I still feel confused. It seems to me that we can directly calculate some or even all the eigenvalues even though the matrix is in relatively huge sizes, for example, $1000\times 1000$, with the formulation $AA^T x=\lambda x$, which can be done in a few second in Matlab. Is there anything wrong with the way I think?

It seems to me that we can directly calculate ... which can be done in a few second in matlab.

Who says that Matlab is calculating it directly, or that it isn't using Gershgorin's circle theorem in its algorithm?  Can you look at the code to see what they're doing?
Finding the eigenvalues of an $n \times n$ matrix is tantamount to factoring an $n^{th}$-degree polynomial.  Factoring any polynomial of degree greater than $n = 4$ is an iterative process, which works better if you can start with a good guess at the roots.
So, when Matlab does an eigenvalue decomposition of your $1000 \times 1000$ matrix, it is not "directly calculating" it -- it's guessing at values, then refining guesses, then reducing the matrix, etc., etc.
Pursuant to that, Gershogorin's circle theorem looks like a handy way to not only do get a set of initial guesses, but to bound the set of possible eigenvalues in a manner that's easy to test for if your algorithm starts to stray.

Is there anything wrong with the way I think?

I would say that believing that typing something into Matlab and getting an answer is the same thing as directly calculating it is an extreme error in thinking.  Someone knows how to solve the problem - you hope.  But if all you know how to do is press buttons in Matlab and get an answer, the set of persons who really know how to solve the problem doesn't include you.


Why are binary signals drawn continuous in time despite being a digital signal The conventional definiton of a Digital signal is as follows:
A digital signal is a signal that is discrete in time and quantized in amplitude. Almost all resources (textbooks, online sources,etc) stick to this defintion and gives emphasis to the point that a digital signal is discrete in time.

A binary signal is a  Digital signal.

But according to the defintion, the signal should be discrete in time to be a digital signal. In the plot, a binary signal is drawn as continuous in time. Then how can this be a digital signal according to the defintion.

A digital signal according to the definition should look like this:

Also I read that a digital signal has different meaning/definition in different contexts.
For example the above mention definition is apt in the signal processing context, whereas in digital electronics where the binary signal mention above is used, a digital signal is a signal that takes only discrete amplitude values(ie. it is quantized in amplitude) and it can be continuous or discrete in time.

I am confused that why is the above definiton given a strict emphasis even in digital electronics context usually. Almost everywhere I look in internet and even some teachers give the conventional definition of a digital signal without mentioning the context.

Are my above findings correct?
I would like to know the exact definitions of a digital signal according to different contexts.

You are right and wrong!

If you would like to treat a binary signal as a discrete-time sequence having values of $0$s and $1$s, then you would plot them as all other discrete-time sequences, using the lollipop display.

However, for those electronic circuits that implement binary logic operators using physical devices (gates etc.), the associated voltage and current waveforms are analog in nature. What's binary there is the information that the voltage or current waveform carries. But the physical existance is analog. So if you want to plot them, then you would do it just like a logic analyser (oscilloscope) would draw them on their phosphorus crt screen...


System identification packages I've been hunting around for an alternative to Matlab's System Identification Toolbox, but the commonality of the terms "system" and "identification" are confounding the search.

Are there any standalone libraries that are written for this purpose?  (in C++, C#, Python, etc.)

The System Identification Toolbox contains facilities for

Linear system identification and determination of transfer function
Non-linear "black-box" models
ODE parameter estimation
Time series identification

One of the more promising substitutes for all of this functionality would be to use R.  Two packages that incorporate some of the time series functionality are the Time Series Utilities and Analysis package and the Flexible Modeling Environment for  Inverse Modelling, Sensitivity, Identifiability, Monte Carlo Analysis.  Also sharing these functions was the Time Series Analysis and Control package

An option for the transfer function calculation and manipulation (which attempts to be a substitute for the Matlab Control Systems toolbox, which overlaps with the System ID one significantly) is the Control Systems Library for Python, but the CSL doesn't have any functions for dealing with system identification, just the ability to manipulate transfer functions (thanks to JasonS for pointing this out).  The Ch Control System Toolkit allows one to model linear systems and determine the transfer function.

I'm sure there are other examples, but they may take some digging to uncover.


1-D Fourier Transform Of A 2-D Image But At An Arbitrary Orientation If one has a 2-D array and would like to take the 1-D Fourier transform along a direction $\theta$ degrees off the horizontal, is there a better/faster way to do this rather than rotating the image by $-\theta$ degrees, taking the 1-D FT across the columns, and then rotating the result back by $\theta$ degrees?

You can take the 1D IFFT using a pseudo-polar grid. Fig. 2 in this paper shows an example of such grid.

Possible choices for the angle $\theta$ are going to be limited. but the transform will be free from interpolation artifacts introduced by the rotation algorithm.


Drawing conclusions from impulse response of a discrete LTI system I have an impulse response of an LTI system which is
h(n)= n*u(n)-u(n-2) where n=[0,3] (1)

Now regarding the stability and causality of this system, i've drawn the conclusion that the system is causal since the output only relies on past or present values of n and it is stable since n is bounded. In this case n=0,1,2,3. Correct me here if i'm wrong but in case h(n) was u(n) - u(n-2) without knowing anything about n, it would stil be stable right? And if h(n) is the same as the original (1) but again without knowing anything about n it would be unstable correct?
Now i want to analyze this h(n) in odd and even components and also find the linear difference equation of h(n) ,but i can't find anything in my textbook that is of any help. It's the first time i've attended a signals and systems course and i'm really struggling to understand a few things, so any help here would be greatly appreciated.

You're right about your conclusions. And now a few hints to help you solve the questions:

Your textbook should state somewhere the definition of even and odd parts of a (real-valued) sequence: $$h_e[n]=\frac12\big(h[n]+h[-n]\big)\\h_o[n]=\frac12\big(h[n]-h[-n]\big)$$

In this case (FIR filter), the difference equation is simply given by the convolution of the input sequence $x[n]$ with the impulse response $h[n]$: $$y[n]=\sum_kx[k]h[n-k]$$


What algorithms can automatically determine a 3D scene from one or a few 2D images? For a project I'm doing, I'm trying to model a scene based on a phone camera's video input, and then insert arrows into it that indicate direction in three dimensions.

I've looked around, and this paper looks quite interesting, and I suppose I could try to adapt its methods:
http://kevinkarsch.com/publications/sa11-lowres.pdf

However, it has the problem that users must input some reasonable spatial constraints into it beforehand, whereas doing this for every frame of a live and constantly changing video is unreasonable.

Another thing that I found allows reconstruction of a 3D object from a few photos of it in various directions, but I can't seem to find the original source anymore.

Any ideas?

SLAM(Simultaneous Localization and Mapping) algorithms can be used to for 3D reconstruction. They offer solutions for both monocular as well as stereo cameras. With single camera they estimate depth with few images and reconstruct the scene. You can find some of the open source solutions here.

Real time 3D reconstruction can done using ORBSLAM and it is also opensource.Here is the related paper for monocular SLAM.


Relation between continuous time transfer function and sampled approximation Suppose I have some continuous time system and associated transfer function:
$$ y(n)=x(n)+x(n-1)$$
$$ H(j \omega) = 1+e^{j \omega (-T_s)} $$
Now suppose I create a discrete-time approximation of this system like this:

The digital system contained within the converter blocks has the transfer function:
$$ H(z) = 1 + z^{-1} $$
I know that the actual transfer function from $V_{in}(t)$ to $V_{out}(t)$ is something else entirely, but is there any meaning in saying that if $H(z)$ is "evaluated" for $z=e^{j \omega T_s}$, then $H(z)=H(j \omega)$? What is the significance of this observation, if any?

$H(jω)$, which is DTFT of the transfer function $h(n)$, is basically a special case of the Z transform $H(z)$, where $z=e^jω$ i.e. $|z|=1$.
In this special case, the unit circle ($|z|=1$) should be included in the region of convergence (ROC) of $H(z)$. In other words, the system (which is LTI in your example) should be BIBO stable.
The system in your example is BIBO stable since it has one pole in $z=0$, therefore its ROC ($z≠0$) includes the unit circle. Hence, we can say that $H(z)=H(jw)$ only when $|z|=1$ i.e. on the unit circle.


why use svd() to invert a matrix? In MATLAB, i compared elapsed time to invert a Hermitian matrix using inverse(), svd(), and chol(). svd() took the longest. So is there any reason to prefer svd() to the other two methods?

The two methods differ, above all, by their applicability to matrix classes.
col (cholesky) decomposes Hermitian, positive-definite rectangular matrices into the product of a lower triangular matrix and its conjugate transpose;
svd (singular value decomposition) factorizes any m×n matrix into the form UΣV*, where U and V are square real or compex unitary matrices, m×m and n×n, respectively, and Σ is an m×n rectangular diagonal matrix with non-negative real numbers on the diagonal.
The method inv internally performs an LU decomposition of the input matrix (or an LDL decomposition if the input matrix is Hermitian), but outputs only the inverse of square matrix only.
Both SVD and Cholesky can be used for computing pseudoinverse of a matrix, provided the matrix satisfies requirement for the method used. The pseudoinverse operation is used to solve linear least squares problems and the other signal processing, image processing, and big data problems.
UPDATE on OP's comment
The matrix can be both Hermitian and not a positive/negative (semi)definite, in which case it is called an Hermitian indefinite matrix. The (generalized) Cholesky decomposes Hermitian positive semidefinite matrices. You can also decompose a negative (semi)definite matrix, say, Anegdef, just call chol()) on -Anegdef, but you cannot compute an inverse of an indefinite matrix with Cholesky, because of unavoidable square roots of negative numbers.
Rank-deficient matrices (having zero singular values) are not invertible, as is the case for most general-case matrices. That is, there exists no matrix A-1 such that A-1A = AA-1 = I, for rank-deficient matrices. Still, many problems which you solve through matrix inversion can be solved for indefinite (and consequently non-invertible) matrices with a generalization of matrix inversion, pseudoinverses of matrices.
I won't bore you with the description of Moore-Penrose pseudoinverse concept, on which subject many articles are readily available online and in textbooks: in short, the matrix A+ is a pseudoinverse of matrix A, if AA+A = A (consequently, A+AA+ = A+). If A has linearly independent columns, the pseudoinverse is a left inverse, because in this case A+A = I, you can check it with a sample low-dimension matrix constructed for this purpose. If A has linearly independent rows, the pseudoinverse is a right inverse, because in this case AA+ = I. In general case, A may have no left/right inverses, but still have pseudoinverse A+ satisfying the equation AA+A = A.
Back to quadratic programming problems: Below, all matrices and vectors are real-valued, so instead of Hermitian conjugation H we have a matrix transpose T here. For an invertible (in a "classic" sense) matrix A, a quadratic form
$$
f(x) = {1\over2}x^{\textrm{T}}Ax-x^{\textrm{T}}b
$$
has a minimum value at xopt = A-1b:
$$
f(x_{opt}) = -{1\over2}b^{\textrm{T}}A^{\textrm{-1}}b
$$
To compute an inverse of A, you use function chol() here.
You can solve a minimization problem for a quadratic form with a non-invertible matrix A, provided A is positive semidefinite, even if A has no inverse in this case. For the singular matrix A, the minimum value is
$$
f(x_{opt}) = -{1\over2}b^{\textrm{T}}A^{\textrm{+}}b
$$
which is reached for all x given by
$$
x = A^{\textrm{+}}b + U^{\textrm{T}}{\left(\begin{align}{0\\z}\end{align}\right)}, \hspace{1cm} z \in \mathbf{R}^{\textrm{n-r}}
$$
, n is the dimension of x, r is the rank of A.
Matrix U in this solution is the factor of singular value decomposition of A = UΣVT, which you compute with the function svd(). So, maybe not all "the intermediate products, whether they be triangular matrices, eigenvectors, diagonal matrices, etc. are irrelevant when the only things that matter are the computation time(FLOPS)"?
In these examples, all matrices and vectors are real-valued, and I am not sure this example can be readily applied to your beamforming problem. Still, SVD is certainly useful in signal processing, e.g., see Image Reconstruction using SVD.
For the sake of completeness: 1) there exists a generalization of Cholesky decomposition useful for computation of a positive semidefinite matrix, and the pseudoinverse can be computed with the help of Greville's algorithm. But this algorithm is implemented through recursion on matrix columns, and the Cholesky speed can be offset by this recursion; 2) the "stability issues" are addressed by the LDL decomposition.
Any way, these algebraic technicalities are not "too much information" for those who work in signal processing: if, even without a preliminary analysis on positive-definiteness of matrices, you start compute using a chol() function and it turns out that you are getting weird results, you now know what might be the culprit.


Time -Domain to Frequency Domain I am constructing a MFCC and come to a section which is kind of confusing me.

It is:

Step 4: Fast Fourier Transform

and has this equation:

Now I understand that FFT needs to be carried out, I just don't understand where the values are meant for:

h(t)

x(t)

H(w)

x(w)

Up to now, I have:

Inputted the raw signal (vector)
Pre-emphasis filtering on the signal
Carried out "Framing" (vector of vectors)
Computed Hamming Window

It's just this step I am unsure of.

Hope someone can help :)

By convention they mean the following:

$x(t)$ = Input in the time domain

$h(t)$ = Filter impulse response (time domain)

$H(w)$ = Filter frequency response

$X(w)$ = Input in the frequency domain

In the equation that you wrote the two "*" symbols mean different things.  In the time domain it is convolution, in the frequency domain it is multiplication.  The equation is just saying that convolution in the time domain (which is how you do time domain filtering) is equivalent to multiplication in the frequency domain.


How to draw the magnitude plot of a variable slope low/high pass shelving filter? (eg. as a function of cascaded shelving filters) I am interested in making a magnitude graph of a variable slope low/high pass filter. I would like to set slope in dB/oct and a cutoff pitch. The y-axis should be amplitude from 0 to 1, and the x-axis should be frequency in Hz.

I have seen people construct very nice variable slope high/low pass filters, where a series of say 16 shelving filters are cascaded to create the desired effect. For example, in Reaktor: https://www.native-instruments.com/forum/threads/3-db-octave-lpf.323734/#post-1611335

These variable slope filters adjust the cutoff frequency of each individual filter to maintain an equally spaced interval between an arbitrary high point and the desired cutoff frequency. They will also evenly adjust the slope of each individual shelving filter to give the final dB/oct slope desired.

Here are the basic magnitude plots of low/high shelves, where c = cutoffFreq, z = zeroFreq, and where I got the equations:

https://www.desmos.com/calculator/ewvksxpd8v
https://db0nus869y26v.cloudfront.net/en/Equalization_(audio)

To graph a variable slope filter in the fashion I desire, I can just multiply ~16 of these variable slope equations in the same way, but I need to be able to figure out the parameters.

To elaborate, I know from the Reaktor variable slope filter project that the cutoff frequency per band in an array of shelving filters would be set from a global cutoffHz/cutoffPitch (where pitch means MIDI 1-135) knob as:

c[i] = cutoffHz * 2^((((135-cutoffPitch)/12) / numberBands) * i)

In an array of 16 bands, i = 0...15 and numberBands = 16.

But I am not sure how to get the z parameters per band to finish the graph. Any ideas for what the z formula might be in this case? Or any other suggestions?

If you can generate a grey noise signal to put through the filter, the amplitude can be plotted as the resulting graph. Also if you pitche a sine input to the function at every frequency that you wish to measure for, the sine's amplitude after the filter will be the result for every frequency.

I understand that i in this graph is the audio sample instance...

*c[i] = cutoffHz * 2^((((135-cutoffPitch)/12) / numberBands) * i)*

Then send sines through it at every of the 16 bands of the EQ and check the resulting maximum, i.e. from a loop of sines of 10 periods of samples at 16 frequencies. perhaps you have to fade the sine amplitude into the filter to avoid an artefact.


Can bandwidth of the closed loop be bigger than bandwidth of individual elements in it? This is a question from control theory but I hope you can help me with it. I was pointed to this community.

So in my control loop I have elements that have cut-off frequency of 3-5kHz and I designed a PID controller for the inner loop and those elements and suddenly I get a closed loop bandwidth of 30-50kHz cut-off frequency.

So the loop is in the figure bellow.

So the PID controller, or voltage controller I designed is a PID with a derivative filter. Current sources are modeled as 1st order filters with 3-5kHz cutoff frequency. How is it possible to get such a high bandwidth of the closed loop if the inner elements have a lower cut-off frequency?

Thanks

Yes, that's perfectly possible and to be expected.

Below is a toy example that starts with a Butterworth low pass filter and puts it in the forward path of a feedback loop. The feedback path has a gain of two.

This makes the whole loop unstable, but it also makes the bandwidth very high.

The original filter's frequency response is below.

And this is the closed loop response.

R Code Below
    #30035

library("signal")

bf <- butter(5, 0.1)
freqz(bf)
dev.copy(png, 'Q30035/Q30035-Butterworth.png')
dev.off()

k <- 2
bf2 <- bf
bf2$a <- bf$a + k*bf$b
bf2$b <- bf$a*k
freqz(bf2)
dev.copy(png, 'Q30035/Q30035-Feedback.png')
dev.off()


Narrow band Interference affecting timing acquisition in spread spectrum communication systems It is often said that a narrow band "in band" interferer will collapse and its BW will spread when it is correlated with a spreading code at the receiver of say a CDMA system while the signal of interest will pop up with a narrower BW (the data BW), and thus could be filtered and that is what is called processing gain.

However, in order for this to take place the code at the receiver needs to be aligned in time with the incoming signal of interest.

How can timing acquisition be implemented in the presence of the interferer? Traditional symbol timing synchronization (if applied to the chips) and which relies on the output of a matched filter would lead to errors in the presence of an inband interferer.

The CDMA Receiver would give you the jamming resistance as you described it. A certain post correlation SNR is required for acquisition per system design requirements (for example that may be 10 dB SNR). The ratio of the spreading bandwidth versus the data bandwidth gives you the processing gain as you described; so if you had 40 dB processing gain (for example), you could withstand a narrow in-band jammer that is 30 dB stronger that the received spread signal, resulting in a post correlation SNR of 10 dB.

This is further explained below in the context of one method for timing acquisition and tracking in a CDMA Spread Spectrum Receiver, specifically the Early - Late gate method as it is fairly straightforward. (Note we are also concerned with Carrier Tracking, but given that is not your question let us assume that the carrier offset is being tracked, so the waveform of interest for this discussion has no carrier offset).

Assume such a baseband waveform with the carrier offset removed as a complex signal sampled at 2 samples per chip as shown in the diagram below:

There are three correlators shown (Early, Late and Prompt) and each performs a multiply, accumulate and dump, and magnitude detection on the correlated output, resulting in a magnitude versus code offset according to the following response:

Acquisition

In acquisition mode, implying we do not know the timing offset (so could be any offset over the duration of the full symbol), the switch shown in the upper diagram is connected to the Offset Search box, which in this case generates a very small offset that is added to the nomimal frequency control word to the Clock NCO (which nominally sets the chip clock to be Chip2x, or twice the chip rate), resulting in a timing offset ramp (the code starting position will sweep slowly relative to the incoming waveform's synchronized position). A threshold detector on the prompt output (not shown) would make a decision (with all the usual probability of false alarm vs probability of detection considerations) that the signal was acquired, seed the Loop Filter output with the current frequency control word and switch the system to Tracking.

[With that context, now your question]

Consider the case of a tone jammer. The correlator output in the diagram above is the signal output after the processing gain. Not shown is the background noise, but you can imagine the horizontal line to be noisy. A typical post processing acquisition threshold is approximately +10 dB SNR (this is again specific to probability of false alarm vs probability of detection design criteria and can be made much more sensitive by having many successive success detections before we decide to switch to tracking).  If we have a spreading code with 40 dB processing gain (this would required the number of chips to be 10,000), then this would imply that we could still acquire in the presence of a single tone jammer that was in the middle of our band and 30 dB stronger than our signal, as received.

Tracking

Once acquired it is important to understand how we can maintain tracking to the desired time offset. To complete the description, once in tracking an error signal is derived from the Early and Late Correlator Outputs, which are each spaced 1/2 a chip from Prompt.

Below shows the results for the Early, Late and Prompt Correlator outputs when there is zero time error.

And below shows the outputs when there is a timing offset such that the Early-Late computation results in negative error. For a simple first order loop, this error can be simply scaled (for loop gain adjust) and added to the Nominal Frequency Control Word. Higher order loops can also be implemented by design in the loop filter.

Note: We assumed the carrier offset was removed. The Complex Prompt Correlator Output can be used as part of the carrier tracking loop using the cross product frequency discriminator based on estimating the frequency offset from as follows:

$$F_{error} \propto I_2Q_1 - I_1Q_2$$

Where $I_1, Q_1$ and $I_2, Q_2$ represent two successive correlation symbol outputs.


Effect of One Sample Latency in Audio Application I am dealing with a customized processor that has a processing latency of one sample (20.8 microseconds) at 48kHz Sampling rate .
What I would like to know is ...does this latency of one sample have any effect on audio application in any way ..
If yes I would like to reproduce this and see for a real time scenario.

One sample latency is pretty much unavoidable so this is fairly normal. Please note that you need to look at the total system level latency which includes DA and AD conversion. Many audio applications use Sigma Delta converters which have much higher latency than one sample, so your processing latency wouldn't make much of a difference.
The impact of latency depends very much on the specific application. One sample is very low and the only applications that would be impacted by this is anything with a feedback loop (active noise cancellation, active excursion control) and any application that has a parallel analog path (effects loop in an analog mixer for example)
It's not a problem for UI responsiveness, audio/video sync, play/sing/talk along, multi-speaker synchronization, etc.
The most latency-sensitive audio application is active noise cancellation and many dedicated ANC processors run at a significantly higher sample rate, just to reduce latency and increase bandwidth of the feedback loop.


Sifting Property of Shifted Impulse In the SE Chemistry forum, someone posted an interesting question on converting a scaled and shifted delta function into Lorentzian by convolution please see "simulating a molecular spectrum". The OP was observing shifts on the x-axis after the convolution of his deltas' with a Lorentzian. This was obviously due to an incorrectly chosen center of the Lorentzian.
In convolution, the horizontal centroids add. So I suggested that a zero centered Lorentzian should be convoluted with the shifted delta and after that there will be no shift on the x-axis i.e., if the delta function is at x=t, after convoluting with a zero centered Lorentzian, there will be a Lorentzian at t - the desired location.
Another user suggested another alternative. His suggestion was that one can also center the Lorentzian at the mean of the x-axis values, $\mu$, and convolute this mean centered Lorentzian with the deltas. There will be no shifts, and indeed this was observed.
The query is:  Centroids add in convolution. It is understandable that when delta was shifted at x=t, and Lorentzian was centered at zero, so the resulting Lorentzian was centered at t+0=t.
Now, why is the other alternative working? This time the delta is at t, and the Lorentzian is centered at $\mu$, the resulting convolution should be at t+$\mu$, but still the deltas do not shift. What is the mathematical basis of getting the same results?
Here is a MATLAB code for a demo.
%Assigning variable names
x=[6200:2:6700-2]'; % x axis
Int=cat(1, zeros(80,1), 1, zeros(169,1)); %
Delta= x(2,1)-x(1,1); % Sampling interval

% Zero centered Lorentzian parameter;
N=[-(length(x/2)):Delta:length(x/2)-Delta]'; % x-axis for the zero centered Lorentzian
A= 1; %amplitude
W= 10; %full width at half maximum
M=0; %centered peak

L1= A./(1+4*((N-M).^2/W.^2)); % Lorentzian equation

% Mean centered Lorentzian parameters
mu=mean(x)
L2=A./(1+4*((x-mu).^2/W.^2)); % Lorentzian equation

F1=conv(Int, L1, 'same');
F2=conv(Int, L2, 'same');

figure (1)
subplot(1,3,1)
stem(x, Int,'.');
xlabel('time'); ylabel ('signal'); title ('shifted impulse')
subplot(1,3,2)
plot(N,L1)
xlabel('time'); ylabel ('signal'); title ('zero centered Lorentzian for convolution')
subplot (1,3,3)
plot(x,F1)
xlabel('time'); ylabel ('signal'); title ('impulse*zero centered Lorentzian')

figure(2)
subplot(1,3,1)
stem(x, Int,'.');
xlabel('time'); ylabel ('signal'); title ('shifted impulse')
subplot (1,3,2)
plot(x, L2)
xlabel('time'); ylabel ('signal'); title ('"mean time" centered Lorentzian for convolution')
subplot (1,3,3)
plot(x,F2)
xlabel('time'); ylabel ('signal'); title ('impulse*mean time centered Lorentzian')

Thanks.

It would be better to show that property using pencil and paper, because you're fooling yourself with that Matlab code. Of course we have
$$\delta[n-n_0]\star g[n]=g[n-n_0]\tag{1}$$
Now define $h[n]=g[n-n_1]$. According to $(1)$ we have
$$\delta[n-n_0]\star h[n]=h[n-n_0]\tag{2}$$
And, consequently, with $h[n]=g[n-n_1]$,
$$\delta[n-n_0]\star g[n-n_1]=g[n-n_0-n_1]\tag{3}$$
No need to waste computing power to show this.


Discussion of simulation results of a signal sequence: amplitude I need your help in understanding if the simulation result is correct.
I have simulated a random BPSK sequence that was upsampled and filtered with a raised cosine filter.
The length of the filter is 150 taps, sampling frequency 16, cut-off frequency 0.5 and transition bandwidth is 0.2.
The noise was simulated with EbN0 = 10 with the following  noise power:
SignalEnergy = (trapz(abs(filtered_signal_tx).^2))*(1/Fs);
Eb = SignalEnergy/(2*Nb);
N0 = Eb./(10.^(EbN0/10));
NoisePower = 2*N0*Fs;

The signal has a phase offset(shift): exp(1j*(2*pi*t+phase_offset));
As a result, I have got a signal, which has a small amplitude. It is less than 0.08...
I don't know if it is wrong or correct?!
What do you think?

Phase offset and time delay (time offset) are NOT the same and I think this may be a primary source of confusion (Given the OP's other recent posts dealing only with time offset but here the formula used introduces a carrier phase offset). Hopefully the following diagrams help distinguish the two.
No Time or Phase Error
Below is the constellation diagram and eye diagram for the (raised cosine) QPSK waveform in the receiver after the matched filter, with no phase or time offset and showing the decision locations in red:

The eye diagram for the real component of the signal is shown below and should be clear that the imaginary component would look identical:
The blue circles are samples with 4 samples per symbol, one of which is at the ideal sampling location (at sample number 3 and 7 specifically; this eye diagram spans over two symbols).

TIME OFFSET
Next a small time offset (timing error) of about 1/3 of a sample is introduced to show the effect of time offset alone.

PHASE OFFSET
With the time offset back to 0, a phase offset of $\pi/10$ is introduced showing it's effect on the constellation and eye diagram.

Conclusion
Phase Offset and Time Offset are two different parameters that each independently need to be corrected (one can exist without the other). Further, phase offset is typically changing with time, which specifically then when changing at a linear rate is a static carrier offset since frequency is the derivative of phase with time, and is corrected with a carrier recovery loop (while time offset is corrected with a timing recovery loop).


How do I retrieve texture using GLCM and classify using SVM Classifier? I'm on a project of liver tumor segmentation and classification. I used Region Growing and FCM for liver and tumor segmentation respectively. Then, I used Gray Level Co-occurence matrix for texture feature extraction. I have to use Support Vector Machine for Classification. But I don't know how to normalize the feature vectors so that I can give it as an input to the SVM. Can anyone tell how to program it in Matlab?

To the GLCM program, I gave the tumor segmented image as input. Was I correct? If so, I think, then, my output will also be correct.

My glcm coding, as far as I have tried is,

I = imread('fzliver3.jpg');
GLCM = graycomatrix(I,'Offset',[2 0;0 2]);
stats = graycoprops(GLCM,'all')
t1= struct2array(stats)

I2 = imread('fzliver4.jpg');
GLCM2 = graycomatrix(I2,'Offset',[2 0;0 2]);
stats2 = graycoprops(GLCM2,'all')
t2= struct2array(stats2)

I3 = imread('fzliver5.jpg');
GLCM3 = graycomatrix(I3,'Offset',[2 0;0 2]);
stats3 = graycoprops(GLCM3,'all')
t3= struct2array(stats3)

t=[t1,t2,t3]
xmin = min(t); xmax = max(t);
scale = xmax-xmin;
tf=(x-xmin)/scale

Was this a correct implementation? Also, I get an error at the last line.

My output is:

stats =

   Contrast: [0.0510 0.0503]
Correlation: [0.9513 0.9519]
     Energy: [0.8988 0.8988]
Homogeneity: [0.9930 0.9935]
t1 =

Columns 1 through 6

0.0510    0.0503    0.9513    0.9519    0.8988    0.8988
Columns 7 through 8

0.9930    0.9935
stats2 =

   Contrast: [0.0345 0.0339]
Correlation: [0.8223 0.8255]
     Energy: [0.9616 0.9617]
Homogeneity: [0.9957 0.9957]
t2 =

Columns 1 through 6

0.0345    0.0339    0.8223    0.8255    0.9616    0.9617
Columns 7 through 8

0.9957    0.9957
stats3 =

   Contrast: [0.0230 0.0246]
Correlation: [0.7450 0.7270]
     Energy: [0.9815 0.9813]
Homogeneity: [0.9971 0.9970]
t3 =

Columns 1 through 6

0.0230    0.0246    0.7450    0.7270    0.9815    0.9813
Columns 7 through 8

0.9971    0.9970

t =

Columns 1 through 6

0.0510    0.0503    0.9513    0.9519    0.8988    0.8988

Columns 7 through 12

0.9930    0.9935    0.0345    0.0339    0.8223    0.8255

Columns 13 through 18

0.9616    0.9617    0.9957    0.9957    0.0230    0.0246

Columns 19 through 24

0.7450    0.7270    0.9815    0.9813    0.9971    0.9970

??? Error using ==> minus
    Matrix dimensions must agree.

The input images are:

Are you using Matlab?  If so then you would either need the Bioinformatics Toolbox, that includes an SVM classifier, or you can download libsvm, which has Matlab wrappers for training and testing.

Then you are going to need some labeled data. Are you classifying liver tummors as opposed to healthy liver?  Then you would need images of liver tumors and healthy livers, each labeled as such.

Then you need to compute some features. What those are, depends on the nature of the problem.  Texture features seem like a good start.  Consider using co-occurrence matrices or local binary patterns.

Edit:
As of the R2014a release there is a fitcsvm function in the Statistics and Machine Learning Toolbox for training a binary SVM classifier. There is also fitcecoc for training a multi-class SVM.


Computing main lobe width of Hamming window in Matlab by plotting I wanted to verify the fact that main lobe width of a Hann/Hamming window is twice that of a rectangular window, using Matlab. I computed the analytical expressions for the respective frequency responses, and plotted them. I do not get the expected result. Can someone explain this?

Code :

omega = -pi:0.01:pi;
N = 20;
fn = exp(-j*omega*(N/2 - 0.5)).*((sin(omega*N/2))./(sin(omega/2)));
shift = (length(omega)-1)/(N-1); % 2pi/(N-1); length(omega) = 2*pi
fn2 = 0.54*fn' + 0.23*circshift(fn',ceil(shift)) + 0.23*circshift(fn',-ceil(shift));
plot(abs(fn)/max(abs(fn)))
hold on
plot(abs(fn2)/max(abs(fn2)),'r')

This is the plot :

As you can see, neither is the main lobe wider, nor are the side lobes smaller.

The problem is the way you shift the transform of the rectangular window to compute the transform of the Hamming window. You can't just shift on the discrete frequency grid, but you have to evaluate the (continuously) shifted functions on the grid. If $W_r(\omega)$ is the transform of the rectangular window

$$W_r(\omega)=\frac{1}{N}e^{-j\omega (N-1)/2}\frac{\sin(N\omega /2)}{\sin(\omega/2)}\tag{1}$$

then the transform of the Hamming window can be written as

$$W_h(\omega)=\alpha W(\omega)-\frac{\beta}{2}\left[W_r(\omega-\omega_0)+W_r(\omega+\omega_0)\right]$$

with $\alpha=0.54$, $\beta=0.46$, and $\omega_0=2\pi/(N-1)$.

In Matlab/Octave you get:

w = -pi:0.01:pi;
N = 20;
w0 = 2*pi/(N-1);
Wr = 1/N*exp(-j*w*(N-1)/2).*sin(N*w/2)./sin(w/2);
W1 = 1/N*exp(-j*(w-w0)*(N-1)/2).*sin(N*(w-w0)/2)./sin((w-w0)/2);
W2 = 1/N*exp(-j*(w+w0)*(N-1)/2).*sin(N*(w+w0)/2)./sin((w+w0)/2);
Wh = .54*Wr -.23*(W1+W2);

plot(w/pi,20*log10(abs(Wr)),w/pi,20*log10(abs(Wh)/max(abs(Wh))))
axis([-1,1,-60,0]), grid

EDIT:

Another error I've mentioned in the comments is the sign of the shifted windows in your formula. The positive sign that you use is valid for a zero-phase Hamming window, i.e. a window that is centered at $n=0$. However, your rectangular window starts at $n=0$, so it's not a zero-phase window. Consequently, combining the two windows as if they were both zero-phase doesn't work. That's where the negative sign in my formula comes from (both windows start at $n=0$). If you change the sign in your code the result is much closer to a Hamming window than before, but it's still off due to the shift on the grid I've mentioned above. In the figure below you can see the difference (blue: exact Hamming; red: your solution with corrected sign):


Implementation of IIR Filter in dsPIC 33EP - Scaling I have to implement an IIR controller in dsPIC 33EP (16-bit, fixed-point, two's complement, wrap-around, 40-bit accumulator, 32-bit multiply, 12-bit ADC) and I am using cascaded, Direct Form I second order sections (ordering up).

Now, I am having trouble with the scaling procedure, since I'read a lot about $L^2$-norm, but can't quite figure out how to apply it properly.

When I get the sos matrix and gain using MATLAB (by either using tf2sos, fdatool, or filterbuilder) using $L^2$ norm, the gain values are always 1 for all of the intermediate sections.

Is this correct?
Or do I need to scale the input for each section?

(I have no experience implementing IIR systems and your help would mean a lot.)

Ben

The tf2sos function takes an input filter of order
$N$, given by $H(z)=\frac{\sum_i^{N} b_iz^{-i}}{\sum_i^{N} a_iz^{-i}}$ and returns coefficients for $N/2$ second-order filters $H_k(z)=\frac{b_{k0}+b_{k1}z^{-1}+b_{k2}z^{-2}}{a_{k0}+a_{k1}z^{-1}+a_{k2}z^{-2}}$ and a gain $g$ such that

$$
H(z) = g \prod_kH_k(z)
$$

So, regardless of what normalization you want from the function (none, inf or two), the above always holds. So to answer your questions:

the gain of each filter section is different (i.e. not necessarily one). (But, you dont care about that)
you do not need to adapt the gain between the filter stages. What comes out of the first stage is directly sent to the second stage. Then, after all stages were processed, you multiply the outcome of the last stage with $g$, what you receive from the tf2sos function.


Question about vanishing moments in wavelet transforms I am reading the book Noise reduction by wavelet thresholding by Maarten Jansen. About vanishing moments, it reads

  To create a really sparse representation, we try to make coefficients
  that live between point s of singularities as small as possible. In
  these intervals of smooth behavior , the signal can be locally well
  approximated by a polynomial.

I don't understand the sentence. If the number of vanishing moments is large, does it mean that the filtered image will be smoother? If I hope to keep edges in images, do I need to use Haar wavelets?

The initial motivation of a sparse representation for a class of signals or images is to pre-separate components with different morphological behaviors.

A traditional model for images is the augmented cartoon: piece-wise smooth + contours + geometrical textures + "noise (unmodelled)".

The mother wavelets (gradient-like operators) can capture isolated singularities at different scale with a few coefficients (and if they oscillate, a bit of textures as well). If the singularities are mixed with a polynomial trend, a standard gradient may also produce coefficients. With a sufficient number of vanishing moments, wavelets somehow "pre-separate" polynomials and singularities, and the polynomial part is cast, subsampled, to a small set of approximmation coefficients.

All in all, you can get few coefficients for singularities, and few coefficients for polynomials between singularities (a long polynomial can be encoded into a few polynomial factors). But still you preserve all the information, especially about the gradients you want to keep.

Wavelets are, under some theoretical conditions, quite optimal for such 1D signals. In practice on 2D data, other aspects come into play and counterbalance the need for vanishing moments: the nature of singularities in images, their finite size, the other wavelet properties, etc.

Resultingly, 1D wavelets work ok for not-too-complex images, but become limited with complicated data, complicated processing, or  high noise conditions.

My experience is that, unless you have specific needs, 2-4 moments suffice in most of the cases, meaning that if finally you use a wavelet with more moments, it will be a by-product of other wavelet properties, and not a strong requirement. There are many options in finding the optimal multiscale  geometric representation.


Improving Main Lobe Width of FFT It is common knowledge that the duration time samples are collected is inversely proportional to the width of the main lobe in the Fourier domain. The simple example is the rectangular pulse in time whose Fourier transform is the sinc function. The width of the rectangular pulse is inversely proportional to the position of the first zeros of the corresponding sinc function.
As a result, if you desire to decrease the width of the main lobe in the frequency domain, you should take additional samples in the time domain.
I am working on a project wherein I can collect only a certain number of sequential samples, let's say $N$. After I collect $N$ time domain samples, I must wait a certain duration before I can collect another set of $N$ samples. Let's define the number of missing samples in the "gap" as $M > N$. As a result, the time samples might look something like the following:
There are 3 scenarios to consider:

"Sample sets": two sets of samples over the time indices $n \in [0, N-1] \cup [M + N, M + N - 1]$
"Aggregated Samples": single sets of samples over the time indices $n \in [0, 2N - 1]$
"Fully Sampled": single sets of samples over the time indices $n \in [0, 2N + M - 1]$

Considering the scenario 1, shown in the first figure, I have collected $2N$ samples but there is a spacing between them. My goal is to recover the Fourier transform of either scenario 2 or 3. In either case, we should see an reduction in the main lobe width in the frequency domain. However, because of the large gap between sample sets in scenario 1, if I use a nonuniform fast Fourier transform (NUFFT computed in MATLAB) or fill the gap with zeros (which is mathematically identical to the NUFFT in this case), the resulting grating lobes dominate the spectrum. Below is a comparison of scenario 3 and scenario 1. As evident in the figure, the frequencies that are closely spaced are considerably distorted due to the finite size as compared with the ideal case. A similar phenomenon can be observed if we were to compare scenarios 2 and 1.

If I were to recover the frequency domain representation of scenario 3, I essentially would be interpolating between the sample sets.
If I were to recover the frequency domain representation of scenario 2, I essentially would be moving the second sample set next to the first and resolving the phase difference between the set of samples that would exist at $n \in [N, 2N -1]$ and the set of samples I have which were taken at $n \in [M + N, M + N - 1]$.
Some other notes:

The signal I am sampling consists of a sum of sinusoidal signals each of which have different frequencies and are likely to be close to each other in frequency, but the number of different components is typically large $>> N$.
For various reasons, I cannot use an autoregressive (AR) method, MUSIC/ESPIRIT/etc. method, Prony's method, or matrix pencil algorithm
I have attempted various compressed sensing (CS) methods with limited success as my samples are taken with such a specific structure rather than somewhat randomly spanning the entire time duration as common in many CS problems.
I have attempted with some success a fully convolutional neural network approach to learn to fit the NUFFT of scenario 1 to the FFT of scenario 2 or 3.

I would like to get the community's thoughts on the problem and possible solutions. Is there anything I am missing that I could use to solve this problem?
My specific question however is this. I believe the problem appears to be somewhat nonlinear and a data-driven approach may be the best option. In this case, what would be an optimal way to exploit the domain specific knowledge of signal and scenario towards choosing the best deep learning approach? What are some options, preferably in the form of research papers, for taking advantage of the scenario using deep learning? A sparse CNN? A generative model (maybe GAN) for interpolating the gap?
Thanks!

Given the phase coherence over non-blanked intervals, a simple and effective approach would be to additionally window the non-blank intervals with a matching window corresponding to the length of the time samples in between blanked intervals. Since the window only tapers the amplitude and does not modify the phase, there will be no issue with different window durations. This will significantly decrease if not eliminate the grating while provide phase coherent processing gain in the DFT from interval to interval over all usable samples. The Kaiser window would be a good window choice for this application. Consider how the DFT for any given frequency bin, as a correlation, de-rotates a signal that may be at that particular frequency such that each sample is aligned in phase, and then the result is summed; our best strategy would be to maximize the information from the samples that are there alone, rather than make our own attempts as to what we think is there to fill in the blanks, to then further convince our assumptions (unless we had other additional external knowledge, in which case such an approach could make sense).
The windowing approach would result with multiple various sized windows over the complete data set corresponding to intervals of non-zero data, and then an FFT would be performed on the entire data set.
To optimize this process further, Sparse FFT algorithms could be employed to take advantage of the gaps in time, but the windowing approach to modify the non-zero time domain data first would be functionally similar.


Power Line interference with different Sampling frequency I've been playing around with EEG data from an online database.

If the power line interference (PLI) is at 50Hz and sampling frequency is 64Hz, then according to the Nyquist theorem, the 50Hz PLI should be aliased with 14Hz.

But according to the database header file, the signal from EEG channel is processed by some prefilter, which is a bandpass filter with passband 0.5Hz - 32Hz.

So in the end, the raw EEG data I get from the database, is analog signal sampled at 64Hz and bandpass filtered in 0.5Hz - 32Hz range.

So my question is, do I have to design a notch filter to remove the PLI at 50Hz if the signal has already been bandpass filtered (32 Hz above removed.)

That depends if the bandpass filter was applied before or after sampling. If it was before sampling (i.e. in analog domain) it worked as an anti-aliasing filter and you dont need to filter out the PL. If it was applied in digital domain, you'd need to suppress it, because then the aliasing would occur. Have a look at the FFT of your signal to see if there is a peak or not at 14Hz, which could come from the PL.

(I guess the filter was applied in analog domain, because it does not make much sense to do it after sampling, since the high cutoff is already the Nyquist frequency.)


Sparse recovery, Restricted Isometry Property for ILL-POSED problems if $\mathbf x$ is $N\times 1$ sparse vector, and $\mathbf A$ is an $M\times N$ matrix with $M<<N$, and we measure $\mathbf y=\mathbf{Ax}$, then compressed sensing theory tells us that we can exactly recover $\mathbf x$ from $\mathbf y$ under some conditions for  $\mathbf A$ like RIP.

But what if $\mathbf A$ is an $M\times N$ matrix with $M>>N$ or $M=N$. Now the equation is not under-determined. But the equation can not be solved directly because $\mathbf A$ is now ill-conditioned matrix. To approximate to solution we can apply regularization. Since $\mathbf x$ still sparse we apply $l_1$ reguarization scheme as $$ ||y-Ax||_2^2 + \lambda ||x||_1$$ and find a sparse solution for x where sparsity level is dependent on $\lambda$.

The questions are

1) Do we still need to consider RIP for $\mathbf A$? What effects the acuuracy of the solution of $l_1$ regularization scheme except $\lambda$?

2) Why do we need RIP in compressive sensing and does it have any relation with the values of $M$ and$N$?

3) If $x$ is not sparse does this regularization fail for sure? Or it can still be work like $l_2$ regularization scheme. Maybe better?

The purpose of compression sensing is to sample a signal below Nyquist given that the signal is sparse in a certain transformation. So it only makes sense when  $M < N$. Othwersie if $M>N$ or $M=N$ then what you have is conventional sampling anyways. An important property of RIP matrices is that it is guaranteed to only change the length of any vector  "very little" as long as the vector  is at least $K$-sparse (has at most $K$ non-zero coefficients. By definition its is the following: $$(1-\delta_S)||x||_2^2 \le ||A x||_2^2 \le (1+\delta_S)||x||_2^2$$ where $x$ is the sparse vector and $\delta_s$ is an arbitrary constant close to zero. So yes, $M$ should be less than $N$ for compressive sensing to make sense and infact all the known matrices that follow RIP are known to exist when $M = O(K(log(N/K))$ strictly. Finding these matrices is not trivial and is a NP Hard problem.

So you do not need to consider RIP for A when $M>=N$, because we are now not talking in compressive sensing terms. The $l_2$ terms should now be the objective because convergence can be achieved in the $l_2$ convergence.

If $M>=N$ and $x$ is not sparse then the regularization could still work but $\lambda$ should be adjusted accordingly with apriori knowledge of $x$ not being sparse. Although I dont see any point of the $l_1$ being part of the objective.


Where can I get the file l1eq_pd.m? To recover the original signal when compressive sensing is used, we get an error, to avoid the error, we need to modify the l1eq_pd.m file. How to go about it? Where can I get this file in Matlab?

First hit on Google seems to be a reasonable start: https://github.com/gpeyre/matlab-toolboxes/blob/master/toolbox_sparsity/toolbox/l1eq_pd.m


Help in designing filter for DUC (Digital Up-Converter) and DDC (Digital Down-Converter) using the built-in DUC and DDC blocks in MATLAB Simulink Can anyone please help me in simulating a simple BPSK transceiver system that correctly receives a signal through a DUC and DDC? I have attached the screenshot of my Simulink system along with constellation diagrams, spectrums and DUC and DDC filter responses. Below are the specs used in the system:

Binary Generator Specs

Probability of zero - 0.5

Initial seed - 61

Sample Time - 1/10e3 (10 KHz Sampling Frequency)

Samples per frame - 44

Digital Up-Converter Specs

Interpolation Factor - [1 16 2] (Total Interpolation factor - 32)

Here, 3x1 Matrix elements indicate the interpolation factor of the three individual stages of interpolation. Simulink DUC block uses a FIR interpolation filter, followed by a CIC filter and a CIC compensator.

Minimum order filter design - Yes

Two-sided Bandwidth of the input signal - 2e3 (2 KHz)

Source of stopband frequency - Auto

Passband ripple of cascaded response - 0.1 dB

Stopband attenuation of cascaded response - 80 dB

Type of oscillator - Sine wave

Center Frequency of output signal - 100e3 (100 KHz)

Input sample rate - 10e3 (10 KHz).

Digital Down Converter Specs

Decimation factor - [2 16 1] (Total decimation factor - 32)

Minimum order filter design - Yes

Two sided Bandwidth of the input signal - 2e3 (2KHz)

Soure of stopband frequency - Auto

Type of Oscillator - Sine wave

Center Frequency - 100e3 (100 KHz)

Input Sample Rate - 10e3*32=320 KHz

As you can see, I do not get correct BER, which means my signal is getting aliased as I receive through the DDC. As seen from the spectrums, I receive correct frequency translation and correct output sampling rates, which means the signal is getting correctly up-converted.

Please help me in debugging where the problem lies. I presume the problem is with the filtered response. The simulink block allows to modify the filter orders of the FIR interpolator, CIC interpolator and CIC compensator. If anyone could help me in correctly design the system, I would be beholden to you. Thank you

Help design my system would be far more detail than what we typically provide here, but I can give you some suggestions that may help or at least lead to a more concise question.

  I do not get correct BER, which means my signal is getting aliased as
  I receive through the DDC.

Not getting the correct BER does not necessarily conclude this.

In your case the error rate appears to be 50% which suggests your receive data is completely independent of your transmitted data - one possible reason this is if the delay for all the filtering hasn't been properly accounted for when comparing the transmit to the receive signals. Instead of comparing the constellations, consider comparing the baseband modulated signals in the time domain as waveform amplitude versus time to help determine what the actual failure mechanism is. That will show you if everything is fine and you just have the expected time offset, or otherwise give you further visibility into where exactly the problem is occurring (look at the in/out of each block in the time domain).

Once you isolate the block that is not performing as you would expect it to, it will either be quite clear to you what the problem is or you will have a more concise question to ask here.

Good luck!


comparison between Fourier transform and autocorrelation function If the Fourier Transform of a signal is enough to determine its spectral density, why do we complicate our lives by defining the autocorrelation function?

Signals can be classified in many ways, and one of them is according to their nature such as being deterministic or random (stochastic). When a signal is deterministic, its spectral content is given by the Fourier transform provided that the signal is absolutely (or square o.w.) integrable; i.e., its Fourier transform exists.

On the other hand when an signal is stochastic, then ideally by definition its Fourier transform does not exist. (a random signal will always have infinite energy as it cannot decay to zero when time goes to infinity because then it will not be random otherwise.)

Therefore for random signals, the Fourier transform will not exist and their spectral density cannot be computed that way. Nevertheless, for such signals an indirect description, based on an autocorrelation measure, can also be defined, whose advantage being that the Fourier transform of it may exist. Therefore the PSD of a random process is defined as the Fourier transform of its ACF.


FFT giving a huge magnitude of first frequency and pretty much zero after that So I have a (visually) very noisy time series signal and I have applied the fast Fourier transform using numpy's fft function. I am wondering why I am seeing the magnitude of the coefficient for the first frequency very large, actually equal to the number of samples ( 256 samples ) in the data, and the rest of the coefficients very low.

I am assuming this means there there is no periodic signal being extracted from the data, but I wanted to get a second opinion about it. How do I interpret this result? Maybe I am doing something wrong.

Input Signal:

Frequency Coefficients:

Zoomed Picture:

That's correct. It's simply a property of your signal: your signal has very small "wiggles" with a very large mean value. Your mean value is about 1 and your wiggles are around 0.01 we would expect the first frequency to be about 100 times larger then all the frequencies.

The first frequency represents the DC component (at 0 Hz), it's simply the sum of the signal.


Best way to do accurate amplitude measurements using FFT To do accurate amplitude measurements on an digital radar receiver, I am searching for the best method to get the amplitude out of a FFT. The signal is windowed using a Blackman window.
So far I'm only looking for the best way to reduce scalloping loss. I'm already in the frequency domain and want to find the real peak of the FFT spectrum. I already  tried parabolic interpolation, but this gives errors if the signal is near the noise level.
What is the best way to reduce scalloping loss errors without using a flattop window?

You can interpolate using the transform of your chosen window function as the interpolation kernel (similar to Sinc reconstruction), and then use successive approximation to find the peak of the interpolated spectrum.
Pre-constructing a polyphase table using an extremely oversampled transform of your chosen window might or might not assist in this computation.


Why the spectral coherence is unity for all frequencies between single-frequency time series and itself In the example below, I am plotting the coherence between time series and itself. The time series do has one frequency.The coherence magnitude was one for all frequencies. I wonder why it is not zero at all frequencies except at the time series frequency. I found that if I include random variable, I get zero coherence at most frequencies except the time series and its surrounding frequencies.
Here is my questions:

Why the coherence is unit for all frequencies however the time series have only single frequency.
Why I got more realistic (coherence values at the time series frequency, and zero else where) values when I add noise to the time series.
Why the coherence (with the noise added case) show more unity for more that one frequency. It looks like that the coherence have unit at the time series frequency and the surrounding frequencies.

import matplotlib.pyplot as plt
from scipy import signal

nt=10000
w1=0.02

xp_pre=np.zeros((nt))
yp_pre=np.zeros((nt))

for t in np.arange(nt):
    xp_pre[t]=np.cos(2*np.pi*w1*t)
    yp_pre[t]=np.cos(2*np.pi*w1*t)

#xp_pre=xp_pre+np.random.randn(nt)*0.03
#yp_pre=yp_pre+np.random.randn(nt)*0.03

plt.figure();plt.plot(xp_pre);plt.plot(yp_pre);plt.show()

#f,coh = signal.coherence(xp_pre,yp_pre,window=signal.get_window('hanning',200),noverlap=100)
f,coh = signal.coherence(xp_pre,yp_pre)
plt.figure();plt.plot(f,coh,'*')
plt.show()

Thanks

According to the documentation, the coherence between two signals $x(t)$ and $y(t)$ is defined as $C_{xy}(w) = \frac{|P_{xy}(w)|^2}{P_x(w) P_y(w)}$, where $P_x(w)$ and $P_y(w)$ are power spectral density estimates of $x$ and $y$ and $P_{xy}(w)$ is an estimate of the cross-spectral density.

In your first experiment, you put in two identical signals $x(t) = y(t) = \cos(2\pi w_1 t)$. Since $x(t)=y(t)$, what you should expect is that $P_x (w)= P_y(w) = P_{xy}(w)$ for all $w$. Therefore, $C_{xy} = \frac{|P_{xy}(w)|^2}{P_x(w) P_y(w)}= 1$ for all $w$.
You might want to argue that the signal(s) only contains a single frequency. In this case we would expect that $P_x(w) =P_y(w) = P_{xy}(w) =0$ for $w \neq w_1$. However:

This would mean that $C_{xy}(w) = \frac{0}{0} = $undefined for $w\neq w_1$.
This wouldn't happen numerically, as due to rounding errors you'd have some epsilons left that could still make this fraction equal to one
You're not dividing power densities but (Welch) estimates of the power densities of a finite window of the signal, hence $P(w)$ will not be exactly zero anywhere. It will contain estimation artifacts and the truncation effects. Since $x(t)$ and $y(t)$ are exactly equal, the artifacts are exactly equal and the coherence is one.

When you add noise to $x(t)$ and $y(t)$ in your second experiment you're adding "signal content" to the frequencies outside $w=w_1$. Since the noise contributions are drawn independently, they are uncorrelated which means that $P_{xy}(w)$ will indeed be zero there, while $P_x(w)$ as well as $P_y(w)$ will not be. Hence the coherence does go to zero.


Echo Removal from PCM 16 bit arrays (time frames)? I am keen to know the basics to remove echo knowing FarEnd Signal and NearEnd Signal. I have read research papers in which its commonly said:
y(n) = (hT)x
where h is the room impulse response and T is transpose and echo signal estimated is y(n).
Now, what I really want to understand is the computation of room impulse response. Given two-time frame arrays of PCM 16 bit. Let's say:
Near_end_sig = [x1,x2,x3,x4,....,xN] (one time frame)
Far_end_sig = [y1,y2,y3,....,yN] (one time frame)
Here I want to know in this case working on digital acoustic signals (PCM-16 bit arrays of size let's assume 1024). How do I compute room impulse response h and how do I estimate echo signal? My main concern is the computation of room impulse response and the echo signal estimation based on the time frames.
Regards,
Khubaib Ahmad

The canonical way to remove an echo (it can be either electrical echo, or acoustic echo) from a signal is by using adaptive filters. Such an adaptive filter should be able to hold a model of an echo path in the form of an impulse response of this path. Having this model, we can use an original signal to create a copy of its echo, so it can be then subtracted from the signal. Here is an echo canceller:

Lets assume that a communication system has two signals in two directions: $x$ and $s$. Lets also assume that signal $x$ creates an echo $z$. As a result the system gets a signal $s$ mixed with $z$: $s+z$. An echo canceller wants to remove the echo, so it uses an adaptive echo model to create a copy of the echo - $\hat z$ and subtracts it from $s+z$ signal. Since the model can't be precise, instead of $s$ we get $s+e$, where $e$ is a residual echo.
In the most simple case an adaptive echo model is a FIR filter, which adapts during a period, when there is no signal $s$ and this becomes a simple system identification problem: it should match $x$ and $z$ so that $z-\hat z=0$. After adaptation we get $H$ - a vector of filter coefficients, which represent a linear part of the echo path's impulse response. Then during regular communication the filter doesn't adapt and we get $s+e$ like this:
$$
s + e = s + z - H* X ,
$$
where $X$ - is an array of past samples of $x$.
The positioning of the filter depends on a system configuration: it can be at the near end or at the far end. The complexity of the filter strongly depends on characteristics of the echo path: long echo path demands a longer filter, also non-linearity of the echo path determines the required filter capabilities (so that a filter not only includes a FIR filter, but also has IIR taps, implements Volterra series or does something else).
In case of an acoustic echo an echo path is long, nonlinear and changing in time. It is estimated as a room echo response. The filter should preferably be placed at the near end, should include nonlinear component, and should be able to adapt during a conversation (i.e. it has to somehow tolerate the double talk). There are multiple techniques for such things.
As for the specifics of your particular signal, it shouldn't change the basic principle of the echo canceller, since it only changes how the signals $x$ and $s$ are encoded.
Here are several links about echo cancellation from the Vocal Technologies site. I've found their site very helpful, having a lot of useful explanations regarding various aspects of echo cancellation.

Non-Linear AEC
Handling of Nonlinearities in Echo Cancellation
Non-Linear Processing in Echo Cancellation
Double-Talk Detection


More Simultaneous Dynamic Range with Fixed Point or Floating Point? What data path would have more simultaneous dynamic range: 32 bit fixed-point or 32 bit floating-point?  Why?
To clarify "simultaneous dynamic range": This is the ability to discern strong and weak signals at the same time, while dynamic range is the total range between between the weakest signal (sensitivity to still meet a performance metric) and strongest signal (saturation or some other decrease in a performance metric such as 1 dB compression, clipping, etc). One example where we would be concerned with simultaneous dynamic range is in an FFT where we have both weak and strong tones simultaneously at two different frequencies and we want to be able to discern both in the same FFT result.
Our vision and hearing are good examples of the difference between simultaneous dynamic range and dynamic range: see https://en.wikipedia.org/wiki/Dynamic_range
Please preface your answer with spoiler notation by typing the following two characters first ">!"

For the same total bit width, a fixed-point representation has a higher instantaneous dynamic range than a floating-point representation.
 For a value to be simultaneously including information about a strong and a weak signal (where the information is being carried by the sum of the strong and weak signal components), the weak signal has to have a value at least as large as the stronger signal's numerical error floor, given by the format's epsilon. In a floating-point representation, every bit allocated to the exponent similarly reduces the number of bits allocated to the significand, and thus increases the format's epsilon. This correspondingly decreases the ratio between the strong and weak signal components that can be included in a single value.
 For example in a standard 32-bit IEEE-754 single-precision floating point representation, 23 bits are allocated to the significand. As a result the smallest number that can be added to a given number $x$ is of the order of $x/2^{23}$ or approximately $8.6\times 10^6$ times smaller. Note that this is consistent with the C programming language's FLT_EPSILON which is approximately $1.19\times 10^{-7}$. In comparison, using a 32-bit fixed-point representation the smallest number that can be added to a large number using near full scale (e.g. 0xffffff00) is of the order of a single unit (e.g. 0x00000001) which is approximately $2^{32} \approx 4.3\times10^9$ times smaller.